Technical Report 1009 



Structural 
Dynamics Model 
of a Cartesian 

Robot 



Alfonso Garcia Reynoso 



MIT Artificial Intelligence Laboratory 



Structural Dynamics Model 
of a Cartesian Robot 

by 
Alfonso Garcia Reynoso 

BSME Instituto Tecnologico de Veracruz 

(1967) 

MSME Instituto Tecnologico y de Estudios Superiores de Monterrey 

(1969) 

Submitted to the Department of 

Mechanical Engineering in 

Partial Fulfillment of the 

Requirements for the Degree of 

Doctor of Science 
In Mechanical Engineering 

at the 

Massachusetts Institute Of Technology 

October 1985 
©Alfonso Garcia Reynoso, 1985 



The author hereby grants to M.I.T. permission to reproduce 
and to distribute copies of this thesis document in whole or in 
part. 



Signature of Author 




Certified by_ 




Department of Mechanical Engineering 

November 2, 1985 



ffjUAu JflMM* 






Accepted by_ 




rofessor Warren P. Seering 
Thesis Supervisor 



/ 



Ain A. Sonin 
Chairman, Departmental Graduate Committee 



Abstract 

Methods are developed for predicting vibration response characteristics 
of systems which change configuration during operation. A cartesian robot, 
an example of such a position-dependent system, served as a test case for 
these methods and was studied in detail. 

The chosen system model was formulated using the technique of Compo- 
nent Mode Synthesis (CMS). The model assumes that the system is slowly 
varying, and connects the carriages to each other and to the robot structure 
at the slowly-varying connection points. The modal data required for each 
component is obtained experimentally in order to get a realistic model. 
The analysis results in prediction of vibrations that are produced by the 
inertia forces as well as gravity and friction forces which arise when the 
robot carriages move with some prescribed motion. 

Computer simulations and experimental determinations are conducted 
in order to calculate the vibrations at the robot end-effector. Comparisons 
are shown to validate the model in two ways: for fixed-configuration the 
mode shapes and natural frequencies are examined, and then for changing 
configuration the residual vibration at the end of the move is evaluated. 

A preliminary study was done on a geometrically nonlinear system 
which also has position-dependency. The system consisted of a flexible 
four-bar linkage with elastic input and output shafts. The behavior of the 
rocker-beam is analyzed for different boundary conditions to show how some 
limiting cases are obtained. A dimensional analysis leads to an evaluation 
of the consequences of dynamic similarity on the resulting vibration. 
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Title: Associate Professor of Mechanical Engineering 
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Chapter 1 



INTRODUCTION 



1.1 Summary 

A position-dependent system in vibration analysis is defined as a system 
whose vibration characteristics are changing as its components change their 
relative positions. Machines and structures such as cranes and bridges are 
examples of such systems. One such system which motivated the present 
analysis, is a cartesian robot (seen in Figure 1.1) as a structure with moving 
carriages) whose dynamic performance needs to be predicted accurately 
during the motion of the carriages. 

The importance of the present analysis comes from the fact that this is a 
precision assembly robot and any residual vibration at the end of the move 
can badly affect assembly operations with small tolerances. So, it would be 
desirable to have a good mathematical model of how the vibrations evolve 
during the motion so that these vibrations could be actively controlled. It 
could also be used as a design model. 

The mathematical model to be developed should be accurate and fast 
to execute when implemented in a computer. The approach followed to 
develop this model was based on Component Mode Synthesis which is a 
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Figure 1.1: Picture of Cartesian Robot 
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method to connect substructures to predict the dynamic behavior of an 
assembly based on the dynamic behavior of each substructure. This method 
(CMS), traditionally applied to connect parts that form a time- invariant 
system can be extended to position-dependent systems under the following 
assumption: that the system is of slowly-varying coefficients; that is, that 
any appreciable change in the system coefficients occurs after many cycles 
of vibration. 

Experimental modal analyses were performed to determine the dynamic 
behavior of each substructure. We used these results applying CMS to de- 
velop the mathematical model. There were difficulties to be overcome. First 
was the problem of parameter identification to accurately scale the eigen- 
vectors of each substructure so that they were orthonormal. This involved 
special curve-fitting techniques applied to transfer functions obtained by 
test. Another important difficulty was to realize that an ill-condition may 
be present when we use several points to connect 2 substructures, so we 
developed a method to detect this problem and to aid in correcting it. 

The model had to be validated by test, and it was done for two different 
aspects: the quasi-static case and the dynamic case. For the quasi-static 
case we fixed the position of the carriages so that the problem became time- 
invariant and then we compared the modal characteristics as obtained by 
testing the complete robot with the modal data calculated by the model. 
The average error in natural frequencies for the first 17 modes was 6.9%. 
The second aspect of the validation involved running the carriages according 
to a given velocity profile at the motor end, and recording the vibration 
at the end effector of the robot; this was compared to the prediction of 
the model when a prescribed motion at the motor end (equal to the one 
obtained in the test) was imposed. These results showed an acceptable 
agreement between the predicted results and the test results. 
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Figure 1.2: Simply supported beam with moving 2-dof subsystem 

1.2 Previous Research in Modeling of Position- 
Dependent Systems 

The problem of studying the vibrations of a position-dependent system in- 
volves dealing with systems of continuous-parameter elements which are 
described by partial differential equations with changing boundary condi- 
tions as the elements change position. This situation, already difficult to 
solve for the case of ideal elements, becomes worse when the elements are 
real, that is, when there are all kinds of non-uniformities present. Then, 
an exact solution is very cumbersome and for practical purposes we should 
use a method that gives a fast and accurate solution. 

One practical approach adopted by civil engineers in the case of the 
bridge problem is explained by Biggs [9], where it is tacitly assumed that 
the system is slowly varying, and they use one point of connection to model 
the contact between the bridge and the vehicle. The analysis simpliQes the 
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single-span bridge structure by taking it as a simply-supported beam; the 
vehicle is taken as a 2-degree of freedom system. To derive the equations 
of motion, assuming one mode of the beam only, the method uses the 
modal equation of motion for the beam A„ +u\A n — .\ f M*) — The term 
F, the force between the beam and the 2-degree of freedom subsystem is 
expressed as a function of the acceleration of the mass in contact with the 
beam (y v ), (see Ggure 1.2) the acceleration of the second mass(z) and the 
spring force (k(z - y v )) in the 2-degree of freedom subsystem. A coupled 
set of differential equations in y c and z results. 

Since the problem in hand, as in many other manipulators, can be stud- 
ied under the assumption of slowly- varying coefficients as is explained later, 
we can extend the idea behind the approach for the bridge problem. We first 
notice that in the method mentioned above (bridge problem), the dynamics 
of the 2-degree of freedom subsystem are introduced to the beam equation 
through the condition of equilibrium of forces at the point of connection. 
This idea is broadly used in methods of connection of two or more sub- 
structures, namely the Impedance method [61] and the Component Mode 
Synthesis method [33], where the starting point is the specification that 
the vibration of any two points of connection should be the same. Then, 
the Impedance method substitutes for the motion, a frequency function 
whereas the Component Mode Synthesis method uses a function of time. 
These methods have several advantages over that of Biggs [9]. They can 
handle several points of connection between parts and they can work in 
3-D in a more convenient way than the previous approach. In addition to 
that, they can use data obtained experimentally, which may lead to a more 
realistic model. 

For the more general case of position-dependent systems, where the co- 
efficients may quickly vjvry during the vibration, the solution may be found 
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if a pattern of vibration is repeated for each cycle of the oscillating system. 
In this case there is extensive literature. A section of this literature is re- 
lated to the vibration of linkages, and for some simple cases the differential 
equations of motion reduce to Hill's equation or to Mathew's equation. One 
example is the research done by Capellen [13]. 

A method for analyzing the complete behavior of industrial robotic ma- 
nipulators with flexible links and including the control systems is presented 
by Sunada and Dubowsky [73] in which they use finite element models to 
analyze the links. The method is made computationally efficient through 
the use of Component Mode Synthesis. 

For the case of the robot where each component remains invariant but 
the boundary conditions are changing with position, the method of Com- 
ponent Mode Synthesis may be more adequate for giving a mathematical 
model which can run fast in the computer since it requires extensive pro- 
cessing of the experimental data to be done beforehand. It uses only the 
modal results for the components in the model. 

1.3 Literature Survey on Component Mode 
Synthesis 

The use of modal characteristics of each component to derive the approxi- 
mate dynamic behavior of the assembled structure may be accomplished in 
different ways as summarized by Hart, Hurty and Collins [34]. These meth- 
ods can be grouped in three categories. The first is Free-free Modal Synthe- 
sis which is based on free-free modes of the components. Hou [33] reviews 
these methods and introduces a new approach which is adequate for han- 
dling experimentally determined mode shapes. This procedure is applied 
in this study to develop the mathematical model of the robot. The second 
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is Component Mode Synthesis which uses boundary-fixed mode shapes and 
boundary displacement functions. This method was pioneered by Hurty 
[36] where natural modes and frequencies of structural systems are deter- 
mined by an energy method, synthesized from admissible mode functions 
selected for the component mode members of the system. The synthesis is 
accomplished by using equations of constraint that follow from conditions 
of force equilibrium and deflection compatibility at the junctions. Other 
authors have introduced variations to the method [34], [16], [3]. All these 
methods are very effective in reaching convergence when only a few modes 
of the components are available. The third category is Component Mode 
Substitution which uses free-free modal displacements and interface loading 
at the attachment points. It was introduced by Benfield and Hruda [4], and 
it has the advantage of improving the accuracy obtained by the free-free 
procedure. 

The accuracy of the model for the free-free case depends on the num- 
ber of modes used. We have found that for this particular application the 
number of modes required to get a good convergence index can be easily 
obtained by test (10 to 16 modes for each component including the rigid 
body modes). 

Other authors have surveyed the existing literature on Component Mode 
Synthesis, Greif and Wu [26] summarize the work done from 1980 to 1983. 
Craig [17] surveys some of the many variations to the basic procedure of 
CMS which have been extensively developed during the past 20 years for 
application in the modal analysis of complex structures. He points out that 
further research on damping synthesis and the use of experimental data to 
either verify or formulate CMS mathematical models should be conducted. 
A technique for modeling a structure using a severely truncated mode 
set is developed by Hintz [32] where in order to improve convergence he uses 
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modes (constraint modes and attachment modes) based upon a static anal- 
ysis of the component structure response due to imposed interface forces 
and displacements. MacNeal [52] uses statically determined deflection in- 
fluence coefficients as well as residual flexibility of the higher, truncated 
modes to improve accuracy. Rubin [67] introduces a method that employs 
an incomplete set of free-boundary normal modes each adjusted to account 
for the contribution of residual(neglected) modes. This method adds resid- 
ual inertia and dissipative effects to the residual flexibility introduced by 
MacNeal. 

On the use of attachment modes in substructure coupling, Craig and 
Chang [18] make an evaluation of the methods that employ attachment 
modes, and they develop a generalized substructure coupling procedure 
which has good convergence characteristics. Hale [27] increases accuracy 
by an iterative procedure for generating improved substructure trial state 
vectors. 

Regarding damping synthesis, there are a few analytical studies. Kana 
and co-workers [43], [44], [45], [46] developed three methods for calculating 
proportional substructure damping matrices from test data. Hasselman 
[30], [31] used diagonal and off-diagonal terms to develop system damping 
matrices that involve various substructure coupling procedures. 

The use of experimental data to develop CMS mathematical models 
and the experimental verification of the simulations has been reported in 
a few papers. Klosterman and Lemon [47] present a building block ap- 
proach to do substructure coupling using either experimentally obtained 
data or modal data obtained by Finite Element analysis. They applied the 
method to complex mechanical machinery consisting of 5 components, and 
used 6-7 elastic normal modes for the most flexible ones. They emphasized 
that the higher modes that were left out in the truncation had a residual 
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flexibility which had a considerable influence on the system frequencies, so 
they included this effect by calculating the residual flexibility subtracting 
from the modal response representation, near the antiresonances, the con- 
tribution of the resonant modes. Then, they added an element with an 
equivalent stiffness derived from the above residual flexibility to the rest of 
the components. Another paper by Klosterman [48] applies this method to 
a complex automotive system. 

Among other works on the use of experimental data we can mention the 
one by Thoren [75] that describes a technique by which orthonormal modal 
vectors computed from dynamic response data, are used to derive mass, 
stiffness and damping matrices for a discrete model of a distributed elastic 
system. The same problem is solved by Berman, Wei and Rao [7], and by 
Berman and Nagy [6]. Ibrahim [40] presents a technique to use a set of 
identified complex modes together with an analytical mathematical model 
of a structure under test to compute improved mass, stiffness and damping 
matrices, this based on the orthogonality condition. Berman [8] extends the 
method of incomplete models and uses it to obtain dynamic equations of 
motion of a helicopter transmission gearbox from shake test data. Goyder 
[25] deals with structures of high modal density and reviews methods of 
mathematical modeling. Ewins and Sainsbury [22] develop techniques for 
obtaining multidirectional mobility data which must be sufficiently accurate 
and complete for vibration analysis of an assembly of connected structures. 

Considerable attention has been given to the topic of convergence of 
CMS and mode selection. Hurty [37], [38], [39] developed a first-order per- 
turbation technique to determine the contribution of truncated modes in 
the natural frequencies of the assembled system as obtained by his method 
of CMS. Morosow and Abbott [55] determined weighting factors to assess 
which component modes should be included. 
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1.4 Organization of the Thesis 

The thesis introduction describes the problem we axe dealing with and the 
approach followed to solve it, as well as the experiments and calculations 
performed. It discusses the existing literature on the topic as well. 

The second chapter, Proposed Quasi-static Model of a Cartesian Robot, 
introduces details about the method of Component Mode Synthesis and 
explains how the mathematical model was developed for the homogeneous 
case, that is, with no forcing terms. It then focuses on predicting modal 
results for a given configuration of the robot. It also explains the difficulties 
encountered with the orthonormalization of the eigenvectors and with the 
ill-conditioning that may be present . 

The following chapter, Proposed Dynamic Model of a Cartesian Robot 
introduces the forcing terms in the model which are due to the inertia forces 
associated with the acceleration of the carriages and to the gravity forces 
acting on a changing system. At this point, the model can predict the 
vibration of any of the test points in the structure. 

Then we show, in the chapter of Experiments and Results, the experi- 
ments that have been done to validate the model and the degree of corre- 
lation obtained with the model used. 

Finally we include three appendices to give more detail about the data 
and results. Appendix A deals with the mode shape data obtained for each 
of the robot components as well as for the assembled robot. Appendix B 
gives detail about the vibrations of the moving robot, in particular about 
the method followed to derive the forcing function. Appendix C shows 
a preliminary analysis done by this author. This work gave knowledge 
and insights on how to deal with problems of systems that change their 
configuration. 



Chapter 2 

PROPOSED QUASI-STATIC 
MODEL OF A CARTESIAN 
ROBOT 



2.1 Application of Component Mode Syn- 
thesis 

The free-free modal synthesis procedure [33] which we have used in the 
analysis and referred to as CMS from now on, can be summarized as follows: 
Suppose that we have 2 components A and B whose natural frequencies 
[w„], [wj,] and mode shapes [0 a ] >[<&>] &ce known. These components can be 
either free in space or fixed, but the interfaces must be free when obtaining 
their mode shapes. Assuming that the eigenvectors are orthonormal (unit 
modal mass) we have two sets of decoupled equations to represent both 
systems: 

{P*} + K 2 ] (p.) = M {/a} 
{p b } + K] { Ph } = [frf {/„} 
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(2.2) 



where the actual vibration would be 

{z„} = [<Aa] {p«} 

{*b} = [4>b] {Pb} 

If we enforce compatibility at the interfaces, then we have for the connec- 
tion points {i a<; } = {xb c } °r [<£aJ{Pa} = [06c] {Pi>} where the subindex c 
indicates that only the interface points are being used. If we apply these 
compatibility conditions, the original generalized coordinates {p a } and {pt,} 
are overconstrained and we have to reduce the order of the system so that 
we end up with a set of linearly independent coordinates {q}. We can 
partition these vectors and the modal matrices as well 

(2.3) 
(2.4) 

[*-] = [[<] m\ ( 2 - 5 ) 

[*.] = [<] (2-6) 

where the superscripts / and D stand for independent and dependent. 
Then, 



' {pi) 




M 




• = . 




{pf } 




[ (rf> J 



{Pa} = 

M = .{ pi } = { <j» } 



WJ [*£] 



rf 



= [<] M 



(2.7) 



Rearranging we have 



W-l-uzrfc] [<]~ l [K] 




(2.8) 



Finally if we include the complete vector {p} we have 



Pi 
\ P? 

Pi 



I 

-m-'fc] [*£\~ l [<] 

/ 



(2.9) 
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or {p} = [T]{g} . This matrix [T] is the transformation matrix that con- 
nects A and B. 

To obtain the equations of motion of the assembled system we first 
combine the two parts as follows: 

.21 



Pa 
Pb 



> + 



w; 



o 

k 2 ] 



[4>aY ' 


(/.) 


o M« _ 


lA J 



(2.10) 



Substituting {p} = [T]{g} and then premultiplying by [T]' we get: 



[Tf m m + in 



w- 







[T] {<?} = pf 



# o 


< 


fa 


o ft 




. f» . 



(2.11) 



These equations are now coupled and after solving the eigenvalue problem 
of this pseudosystem we get {u/} and [ip]. Finally for the assembled system 
we have that w and 4> are given by 



{"} = {«'} 



(2.12) 



and 



M = 



in w 



(2.13) 



[0a] 

o [h} 

Notice that the name free-free for this method does not imply that the 
components have to be tested on a free-free condition (floating), but rather 
that the points of connection of the components should be free during the 
test (as opposed to the fixed-boundary case), while the other points can be 
either fixed or free. 

One important characteristic of this method (CMS) is that it can ei- 
ther use mode shapes for each component as obtained by a finite-element 
technique or the mode shapes can be obtained experimentally. There are 
complicated cases where we may be better off using test results as data for 
CMS in order to get a good model. Another important advantage of CMS 
is that it handles data in matrix form. This allows us to tackle complicated 



CHAPTER 2. PROPOSED QUASI-STATIC MODEL 



23 



Z, measured from neutral position 




Figure 2.1: Two-dof subsystem moving on beam 



problems with many degrees of freedom at the connections and in 3-D with 
ease. This CMS approach works well for time- invariant cases, that is, when 
the points of connection are always the same and each component remains 
invariant. This approach can be extended to the case of position-dependent 
systems for the case of slowly-varying coefficients as shown below. 

Consider the following example: Suppose we have a bridge problem as 
explained by Biggs [9] and illustrated in Figure 2.1. We will here develop the 
equations of motion for this system using CMS methods and then we will 
compare the results with those obtained by Biggs through use of another 
method. The bridge is modeled as a beam and the vehicle as a 2-degree of 
freedom system traveling at constant velocity as shown in Figure 2.1. 

Assume that M Vu slides without friction and always in contact with the 
beam. Assume also that the time Ti = -, required to travel the distance 
L, is much larger than a natural period of oscillation t n : r n >> t n . We 
can conceive two time scales in this problem; a fast time t that follows each 
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cycle of vibration and a slow time r that follows the rigid body motion of 
the carriages. This variable e is defined as the ratio of the slow time r to 
the fast time t : e = j , and for the beam problem its formula is 

«=f| (2-14) 

or in other words e is the ratio of the period of the dominant mode of 
oscillation to twice the travel time. 
Assume r = et , then x — vr — vte . 



2.1.1 Characteristics of the beam 

Consider the first mode only for a simply supported beam: w A , <f>A — 
JJLsmf,f A = 0. 

2.1.2 Characteristics of the 2- degree of freedom sys- 
tem 

The 2-degree of freedom subsystem, free in a 1-D space has the following 
natural modes: w|, = , u/|, = k v (fejJ&O 

C\ Ci 



[4>b] = 



Ci —Cijf- 



M„ 



(2.15) 



where d = VM ^ +Mvu , and C 2 = yf^ (Mvu+Mv ,) 

M va g 



[/b} = 



(M V3 + M Vil ) g 



(2.16) 



2.1.3 Compatibility condition 

At the connection point we have y VA — y VB or 

l<f>A c }{PA} = [4>D r ] 



PC 
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Substituting we get 



2 . -nvte . . 
— -sin— -{Pa} = 
mL L 






Pn x 



(2.17) 



eliminate pu 2 to get 



Pa 
Pd, 
Pd 2 



= fri 



Pa 
Pd 1 



(2.18) 



where 



[T] = 



1 





1 



With this transformation matrix [T] we can form the equations of motion 
of the assembled system in terms of the remaining generalized coordinates 
Pa 



{?} = 



Pb> 



> to get 



Hm+m<i} = {n 



(2.19) 



where 



h = [rnn [*] = in 



{F} = [T]< 



u 2 A 

w| a 



PI. 














(M va + M„ u ) g 
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Substituting the expression for [T] the equations of motion become: 



+ 



i , _2_M ilSL (h / f , ** \gi_2iruli _M, 

symmetric 
^ + ^Afe(M„ u + M t)8 ) 2 sin 2 ^ 
symmetrtc 

gCM^ + M^y^sin^ 





- V^M^+M^j sin ^f 



1 + 



Mm 
M„. J 



9i 

<72 



jfe\/^( M - + M -)sin 



7T«t< 









(2.20) 



To compare this solution with the one provided by Biggs we must determine 
y c , the vibration at the midspan of the beam, and z , the vibration of the 
mass M V3 . These are obtained by 



Vc 

z 



Substituting the expression for [T] and solving for q we get: 



(2.21) 






m 



LML 



\2{M VU + M VS ) 



( M„ u +M v . 



V W. 

L M„. V mi 



MlLU, 17 gjjj XVtt 



Vc 

Z 



(2.22) 



Also, the second time derivative of q is given by: 



IH 



mLM2. 



12 




2(M„„+M„.) 



i/ 



A*. 



ML 



M v . V mi Sm *Y 



+ 



M v . V rnL L COS Z, 





TTVtt ft 



Vc 

z 



+ 



(2.23) 





Substituting these expressions for q and q and premultiplying by 



Vc 

z 



1 ^-^sin^ 



(2.24) 
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we decouple the mass matrix and finally get 



+ €< 



1 + ^A^sin 2 ^ 
M va 





-M„ u (f) 2 sin=f 



%{M va + M vu )sm 




irvr 
L 



Z 



} + 



"a + ^t sin -r- -nrf sin 



Vc 

z 



) +« 



-Aj, sin sst 




2 jrtiU 2k,, • nvtf 

L tnL Sm L 

k 2 



Mm £ COS L 






(2.25) 



These equations are the same as the ones obtained by Biggs when we neglect 
all terms in e . So, the method of CMS gives the right answer when e << 1 
which is the case for systems with slowly-varying coefficients. 



2.2 Modal Representation of the Robot Com- 
ponents 

In order to use CMS, each component of the robot must be represented by 
its rigid body and elastic modes, and these modes must be tridimensional. 
The robot components are seven as can be seen in Figure 2.2. They are 
called STR,X,Y,Z,MX,MY,MZ whose description is given below. The 
cartesian robot is mounted in the Artificial Intelligence Laboratory of MIT 
in building NE43. It was designed and built by a team of students under 
the supervision of Prof. Warren P. Seering. 

STR: The component STR is formed by the structure that holds the car- 
riages and it also includes the dynamics of the floor on which it is 
mounted, that is the 9th floor of Building NE43 of MIT. This struc- 
ture has on its top part two parallel ways which guide the Y-carriage 
when it moves. This structure is formed by 6 pieces joined together by 
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Y-carriagc 




X-carriage 



Z-carriage 



Figure 2.2: Components of cartesian robot 



CHAPTER 2. PROPOSED QUASI-STATIC MODEL 29 

long and thick bolts so that a large torque can be applied to tighten 
them, thus avoiding any creep-friction problem that would introduce 
nonlinearitics in the vibration analysis. 

Component X: The component X is the X-carriage which moves in the 
X-direction; it is supported by the ways of the Y-carriage and it 
supports, in turn, the Z-carriage. When it moves it is acted upon by 
the motor MX through a lead screw. 

Component Y: The component Y is the Y-carriage which moves in the 
Y-direction; it is supported by STR and it supports, in turn, the 
X-carriage. It is acted upon by the motor MY through a lead screw. 

Component Z: The Z-component is the Z-carriage which moves in the 
Z-direction and is supported by the X-carriage. It is acted upon by 
the MZ motor through a rack and pinion transmission. 

The component STR being fixed to the floor does not have rigid body 
modes, but the rest of the components do have, and to calculate them it 
is necessary to know the mass and inertia properties which are shown in 
Table 2.1. 

Before calculating the rigid body modes, it is necessary to specify the 
points on the structure to be measured. These points must include at least 
the connection points of the components, and for the case of the parts that 
have ways that slide, we take as many points along the ways as possible. 

The rigid body modes are six in general, three rigid body translations 
along X,Y, and Z directions and three rigid body rotations about those 3 
axes. It can be shown that these modes produce translations in the x,y 
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and z directions which are given by the following formula: 



<t>y n 






1 











(Z-Z rt ) 


(y r rg ) 
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v/M 





y/T v 





[X X ra ) 


o 
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(V-K„) 


(X-Xc.) 


o 






y/'M 


v7. 


y/I. 





Sin 
<*2n 
S, n 
Sin 
S$ n 
Sen 

(2.26) 
where M is the mass of the component, I x is the moment of inertia about 
its center of gravity(X C9 , Y cg ,Z cg ) with respect to the X-axis and so forth. 
4> Xn represents the mode shape value of a point with coordinates (X, Y, Z) 
along the X-direction for mode n. <f> Vn corresponds to direction Y, and so 
forth. <5 ln to <5 Cn are Kronecker deltas which have the value of zero unless n 
coincides with the first subscript, in which case they take the value of one. 

The above formulas are correct when the principal axes of the compo- 
nent are aligned with the X,Y, and Z axes but if they are not, then one way 
to correct this is by calculating the mode shape values along the principal 
axes first, and then rotating them to coincide with X,Y, and Z. 

The mass, moment of inertia and position of the center of gravity can 
be either calculated or determined by test. For the moments of inertia 
we can swing the component as a pendulum and measure the period of 
oscillation To, and from it we can calculate the moment of inertia by means 
of I cg = ^MgLTQ — ML 2 , where L is the distance from CG to the center 
of rotation. This method seems to be simple to use but it is very sensitive 
to errors when L is short and the position of the center of gravity is not 
precisely determined. Another procedure consists in using flexible strings 
to obtain a rocking motion as well as a translational motion from which 
the moment of inertia can be determined. The way wc proceeded was to 
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do analytical calculations of the moments of inertia. 

With regard to the elastic modes, each of the components, with the 
exception of STR, was tested on a free-free condition, that is, "floating" 
by mounting it on soft foam pads; then the Impact Test was applied to 
a series a test points(the same used for rigid body modes), and the mode 
shapes obtained. These mode shapes can be expressed as arrays where the 
number of rows is the number of degrees of freedom (3 for each point times 
the number of points) and the number of columns is the number of modes. 

<f>ll • -<t>\n 

4>ml • • 4>rrm 

A pictorial representation of these modes (in three dimensions) was ob- 
tained in a Structural Dynamics Analyzer to check continuity and compat- 
ibility of the modes (to find out if there were bad measurements). 

The structure STR was tested differently since it was attached to the 
floor. Here the modes were obtained by a random excitation test using 
an electrodynamic shaker. The test points along the ways were marked at 
1-inch apart to have good resolution, since the Y-carriage runs over them. 

The frequency range for each of the main components varied since each 
one had its own frequency interval and we wished to have a minimum of 5 
elastic modes of each component. Table 2.2 shows the natural frequencies 
determined for each of the main components. The mode shapes can be seen 
in Appendix A . 

The components MX and MY which are motors coupled to lead screws 
were treated in a special way. MY is mounted on STR and its dynamics, 
other than torsion, were already included in the test done on STR. Similarly, 
MX is mounted on the Y-carriage and it was tested already, except for 
torsion. So, all we needed was to calculate the torsional modes of MX and 
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MY. These were as shown in Table 2.3. The component MZ has a short 
shaft with a pinion at the end, and it was modeled as being rigid in torsion. 

The modal analysis method employed by the Structural Dynamics An- 
alyzer involves curve fitting the transfer functions to determine the modal 
parameters of the component. This technique, which calculates eigenvectors 
by assuming the existence of only one mode at a time, is not very accurate, 
that is, it does not give correct absolute values of the mode shapes. The 
relative values it predicts are probably more reliable. 

In order to normalize the elastic mode shapes (for unit modal mass), a 
convenient step when formulating the equations of motion based on exper- 
imental data, we proceeded as explained in the next section. 

2.3 Modal Parameter Identification 

When a modal test is performed on a structure, we obtain the natural fre- 
quencies and some scaled version of the mode shapes. This information 
represents a family of dynamically similar structures, and it is not until we 
determine the right scaling factor to normalize the mode shapes that we 
have really identified the particular structure under test. This normaliza- 
tion factor is difficult to obtain with accuracy, when working experimen- 
tally, because we have to make assumptions about the modal damping in 
the structure when we use transfer function measurements. So, even though 
there are many algorithms to curve fit the measured transfer functions as 
discussed below, the calculated mode shape values (or scaling factors for 
normalization) will depend on the algorithm used, and there is an impor- 
tant variance in the results. This in turn results in an uncertainty in the 
predicted results of CMS. 

There are a number of ways that people have proceeded to normalize 
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the mode shapes based on experimental data. One simple approach used 
in exact analysis of continuous systems, which may work for structures of 
simple geometry (like a uniform beam for example) consists of applying 
the formula for normalization of a continuous beam: C = < i « 

/v ^„ l where C would be the scaling factor to get the correct mode 

shape {<f>} = C{ij}}. This formula can be extended to other simple shapes 



where we integrate over the contour: C « /v ^„ * „ 

Another procedure to find the correct mode shape values is based on 
multi-excitation of the structure [66], in which the structure is excited at 
several points simultaneously in such a way as to excite one single mode 
of vibration; then the vibration is measured at each of the test points and 
that vibration gives the mode shape directly. This method which implies 
having more sophisticated equipment was not used in our laboratory. 

The most common method used is based on a drive-point transfer func- 
tion measurement [77], which in our case was obtained with high frequency 
resolution, and then a curve-fitting technique was applied to it to extract 
from each mode its modal parameters. One of the formulas for the drive- 
point transfer function is : 

where w is the frequency, A and B are parameters appearing in Au 2 which 
represents the contribution of the higher modes, and B which is the con- 
tribution of the lower modes. Each mode has three parameters: w* is the 
natural frequency in rad/sec, & is the modal damping ratio, and 4>kk is 
the true mode shape value at the drive-point. This formula assumes the 
following idealizations: 

• Damping is linear and viscous. 

• Damping is proportional. 



H(u) = Aw 2 + B+Y ?fc*LJ , 2 2? } 
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In order to get these modal parameters, different algorithms that curve- 
fit the experimental data have been prepared, some of them assume there 
is one dominant mode only, and they compensate for the other modes 
assuming these are far apart. In this single mode group of techniques we can 
mention the circle fit to the Nyquist polar plot (real part vs imaginary part 
with frequency as a parameter), the quadrature picking, and the rational 
fraction least squares [11]. 

Another group of techniques make multimodal curve fits to the trans- 
fer function. In general they involve a large number of iterations. They 
may also have problems of convergence since these are nonlinear paramet- 
ric fits (linear in (f>\ k , nonlinear in a/* and &), and an iterative procedure 
is necessary. Among these methods we can mention the method of ratio- 
nal fraction polynomials [11], the Marquardt Algorithm [54], and Global 
Fitting [24]. Brown, Allemang, and Zimmerman [11] describe parameter 
estimation techniques that can be used to determine modal parameters 
from experimentally measured frequency response or unit impulse response 
with particular details on multidegree of freedom techniques. More [54] dis- 
cusses an efficient implementation of the Levenberg-Marquardt algorithm 
and shows that it has very strong convergence properties. 

The way we proceeded, which was easy to use and readily available, 
was to utilize a Simplex technique that minimizes the sum of the squares 
of the residuals (or the differences between the data and the correspond- 
ing calculated values based on succesively improved parameters). Caceci 
and Cacheris [12] present this technique of curve fitting and discuss its 
application, advantages and disadvantages. The Simplex program used for 
curve-fitting admits any function, linear or nonlinear, whose coefficients are 
the parameters to be determined. We used equation 2.27 with the modal 
parameters to be adjusted a;*;, &, and 4>kk with n varying from 5 to 13 
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Figure 2.3: Curve fit of transfer function 



according to the component. We first tried to determine all parameters 
simultaneously, which are 3n, but the process was converging so slowly (if 
at all) that we could not proceed. Then, taking as initial values for these 
parameters the values obtained by a single-mode curve-fitting technique, 
we applied the Simplex method to the determination of three parameters 
at a time (for one mode), but using equation 2.27 with all the modes of 
interest. Here convergence was achieved in a small number of iterations. 
Then we proceeded in the same way with each of the subsequent modes, 
each time updating the corresponding parameters. This procedure tended 
to converge a little faster. One of the results of the curve fits is shown in 
Figure 2.3. 



CHAPTER 2. PROPOSED QUASI-STATIC MODEL 



36 



nut 



load applied by 
housing 




Figure 2.4: Detail of lead screw and nut 

2.4 Physical Connection of the Components 

Before applying CMS we have to define the degrees of freedom of each 
connection, that is, the vectors representing the motion of the points of 
connection. This depends on how the components are physically connected. 
When the carriages are constructed to form the cartesian robot, they 
are connected by means of preloaded cam followers that run on straight 
bars or ways as illustrated in Figure 2.6 . There are a number of these 
connections between each mating pair of components as can be seen in 
Figure 2.5. With regard to the connection of the lead screws, there is a 
special nut mounted on the component to be moved, with lubricated balls 
inside, thus transforming the rotation of the lead screw into a translation 
of the carriage, (see Figure 2.4). 

Since the cam follower acting on the way exerts a normal force (if friction 
is neglected), we can consider as the degree of freedom in that connection 
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Figure 2.5: Schematic representation of physical connections 
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Figure 2.6: Cam followers on ways 



one coordinate or vector perpendicular to the way. Friction is accounted 
for by putting (in selected positions only) one vector(degree of freedom) 
parallel to the way. The nut that connects to the lead screw is represented 
by a vector(degree of freedom) parallel to the lead screw axis. So, the 
degrees of freedom of all these robot connections are defined as schemat- 
ically represented in Figure 2.7. Notice that we use translational degrees 
of freedom, but not rotational ones. Bending moments and torques are 
transmitted between two components by using several degrees of freedom 
in the connection. 

Friction was included as a friction degree of freedom because the pre- 
dicted mode shapes of the synthesized robot(without friction degrees of 
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L • lead screw-nut degree of freedom 
f - friction degree of freedom 
p - rack-pinion degree of freedom 



Figure 2.7: Degrees of freedom of physical connections 
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freedom) showed a relative motion between the parts which was not ob- 
served in the experimental mode shapes. There is also a more fundamental 
reason for including the friction degrees of freedom; a Coulomb friction 
force opposes motion and tries to equilibrate forces until it saturates, that 
is, when its static limit is reached. This saturation makes its behavior 
nonlinear as opposed to a linear spring, for example, which ideally never 
saturates. So, if we assume that the forces never reach saturation, and that 
is because the cam-followers are preloaded against the ways, then we can 
say that the friction effect is like a bolt joint which does not allow slide 
motion during a cycle of vibration. Nevertheless, to avoid redundancies we 
introduced only a selected number of these friction degrees of freedom just 
to prevent some relative motions between components based on experimen- 
tal results. 

2.5 CMS Model of Connections 

After deciding what distribution of degrees of freedom can represent com- 
ponent connections and having the modal representation of those individual 
components, we can use CMS to establish the compatibility conditions at 
the interfaces and end up with a prediction of the dynamic behavior of the 
robot. Figure 2.8 shows a schematic view of the connections illustrated in 
Figure 2.5 by representing each component with a box and each connection 
with an ellipse. The numbers inside the ellipses are the numbers of degrees 
of freedom used, whereas the numbers inside the boxes are the number of 
modes used. 

The compatibility condition is the starting point for CMS. It says that 
the points of connection of two adjacent components have to move together 
in the same way, that is, x Al = x Dl . A slightly different condition is used 
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no. inside box: total no. of component modes 

no. above arrow: no. of generalized coordinates chosen to be dependent 

encircled no.: no. of generalized coordinates chosen to be independent 



Figure 2.8: Box diagram of Component Mode Synthesis 

when we describe the assemblies MX and MY (motor and lead screw). Here, 
the motion of the component that is moved by the lead screw is equal to 
the sum of the motion of the supporting component plus the translational 
motion of the nut caused by rotation of the screw, that is 1,4, = xq, + xn. 
For example, in Figure 2.9 we show the case of MY, STR and Y. Here, if 
Y moves to the right, it may happen that only MY rotates, or only STR 
translates, or both move. But in any case the sum of motions of STR and 
MY must be equal to the motion of the Y-carriage. Something similar 
happens between the Z-carriage and X and MZ. 

The mode shape data of the components can be defined as in table 2.4. 
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Component 


Mass (Kg) 


CU(riim) 


U,(K<J-m 2 ) 


hv V 


'*. 






i = 359 








X-carriage 


17.556 


y = 451 
z= 130 


0.3846 


0.18006 


0.25495 






i = 185 








Y-carriagc 


25.4016 


y = 471 

2 = 74 


0.8060 


2.0194 


2.1081 






x = 368 








Z-carriagc 


12.2458 


y = 542 
z= -160 


0.7404 


0.7404 


0.0234 



Table 2.1: Mass and Inertia properties of Robot Carriages 
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Figxire 2.9: triple connection 
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Mo(l<: No. 


STR 


X-carriago 


Y-carriagc 


Z-carriage 


1 


14.G7 


168 


121 


587 


2 


21.77 


218 


140 


623 


3 


36.87 


431 


150 


772 


4 


40.89 


606 


181 


1115 


5 


51.78 


- 


200 


1471 


6 


62.62 


- 


234 


- 


7 


81.81 


- 


356 


- 


8 


96.95 


- 


393 


- 


9 


128.37 


- 


415 


- 


10 


136.26 


- 


506 


- 


11 


156.82 


- 


- 


- 


12 


176.83 


- 


- 


- 


13 


193.90 


- 


- 


- 



Table 2.2: Natural Frequencies(Hz) of Main Components 



Mode 


/„(Hz),MX 


<{> n (x)MX 


/ n (Hz),MY 


<M*),MY 


1 





259.56 





254.98 


2 


1321 


346cos(2.734£-4u; n :r) 


993 


298 cos (2.734£ - 4w n i) 


3 


3963 


354cos(2.734£-4a; f ,:r) 


2979 


307cos(2.734£-4w n z) 



Table 2.3: Natural Modes in Torsion of MX, MY. 
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Each component has at least one joint (with several degrees of freedom) 
and for each of this joints there is an associated modal matrix for which 
the number of rows is the number of degrees of freedom, and the number 
of columns is the number of modes. Notice that these modal matrices are 
partitioned in such a way that the first submatrix (sufix R) corresponds to 
the independent generalized coordinates, and the second submatrix (sufix 
O) corresponds to the dependent ones. The number of generalized coordi- 
nates for each component can be seen in Figure 2.8. The encircled numbers 
above the top boxes, and the ones below the bottom boxes are the number 
of independent generalized coordinates. 

Based on these compatibility conditions we can establish the following 
matrix equation: Ap D = Bp 1 , where [^{p^} is 

FSTROl -FYALO 

FYA20 



FSTR02 -FYAZO 

-FYAAO 



' (2.28) 



PSTR 

Px 

P*Y 

V l z 

Pmy 

{ Pmx 
(2.29) 

From this equation we get the dependent elements of the generalized coordi- 
nates p D as a function of the independent ones p 1 . Then, the transformation 
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and [BHp 1 } is 
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Component 


Modal matrix 


Size 


Connection (Figure 2.8) 


Fixed structure (STR) 


[FSTRR1 FSTROl] 


8x 13 


a 




[FSTRR2 FSTR02] 


1 x 13 


d 


X-carriagc 


[FXA1R] 


10 x 10 


b 




[FXA2R] 


10 x 10 


c 




[FXA3R] 


1 x 10 


e 




[FXA4R] 


lx 10 


f 


Y-carriage 


[FYA1R FYAIO] 


8x 16 


a 




[FYA2R FYA20] 


10 x 16 


b 




[FYA3R FYA30] 


lx 16 


d 




[FYA4R FYA40] 


lx 16 


e 


Z-carriage 


[FZT1R FZTIO] 


10 x 11 


c 




[FZT2R FZT20] 


1 x 11 


f 


Ball screw-motor, X-direction 


[FMX1R FMX10] 


1X2 


e 


Ball screw-motor, Y-direction 


[FMY1R FMYIO] 


1X2 


d 


Ball screw-motor, Z-direction 


[FMZIO] 


IX 1 


f 



Table 2.4: Modal Matrices of Robot Components. 
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matrix which is implicitly defined by 

= [T\{P') (2-30) 



P D 



can be formed. So, by following the procedure already described in section 
2.1 we can obtain the equations of motion of the cartesian robot, which 
for the homogeneous case have a solution described in terms of the natu- 
ral modes of vibration. These modes were synthesized in the computer on 
the base of the previous analysis and then transmitted to a Structural Dy- 
namics Analyzer in order to take advantage of its display capabilities; thus 
we were able to see a 3-D representation of each mode shape on the ana- 
lyzer screen, the comparison of experimentally determined with synthesized 
mode shapes could be conveniently done. Before showing some of these re- 
sults, we should discuss some issues about convergence and conditioning of 
this CMS process. 

2.6 Convergence of CMS 

As it has been pointed out by several authors, among them Bathe [2], the 
technique of CMS makes use of the Rayleigh-Ritz analysis [2] in the sense 
that it takes as the Ritz-basis vectors the eigenvectors of the individual 
components and then through a linear combination of them, it gives the 
best solution available that corresponds to those basis vectors. This is only 
an approximation since we are not using all the modes of each component 
but rather we have truncated by establishing a finite frequency range of 
study. 

To understand why CMS is considered a Rayleigh-Ritz analysis and un- 
der what conditions this is true, we perform the following analysis. Suppose 
we somehow know the complete mass and stiffness matrices of the assem- 
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bled system formed by the substructures A and B, but suppose that the 
eigenvalue problem is too large to be solved: 

[*][*] - A[M)[*] (2.31) 

We want to reduce the system order from m to s by using the Ritz- 
basis vectors <f> (m x s). In so doing, we are trying to get an approximation 
to $, call it 4>. Put \4>] — [4>\[x\ where [i] are the Ritz coordinates to be 
calculated. 

Then form the Rayleigh quotient 

m = mm (2 . 32) 

which becomes after substitution of $] 



where [K] = [4>\ T [K] [4>] and [M] = [<A] T [M][<£] . 
Minimize the Rayleigh quotient: 

dp{4>) = 2mS' =1 Zj-kjj ~ 2kZU x,m ty = 
dxi m 2 



(2.33) 



(2.34) 



Use p = 4- to get 

E;=i x^ - pmv) = fori =1,2, ..a (2.35) 

This is equivalent to the eigenvalue problem [K][x] = /j[Af][z] of reduced 
order s which we suppose we can solve. This process parallels the CMS 
process as illustrated in Table 2.5. 

From this analysis we conclude the following: 

• The Ritz coordinates [x] are the same as the product [T]^] in CMS. 
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Concept 


Component Mode 
Synthesis 


Rayleigh-Ritz 
Analysis 


Original 
equations 


\Pa] + [»*a]{pa} = {0} 

\pb] + HKpb} = {0} 

only n modes are known 


[K\[4>\ = \[M\[4>\ 

[K\,[M\ are known, 

but not \4>\. 

Problem too 

big to solve 


Choose 

Ritz-basis 

vectors 


[4>] = MM 


[<t>] , a guess 


Apply 

physical 

principles 


compatibility: 

[4>a c ]{pa} = [4>Bc]{Pb} 

Eliminate overconstraining: 

{p} ->{<?} through {p} = [T]{q} 


assume [0] = [</>][z], 
minimize p(<f>), 

or minimize 
potential energy. 


Reduced 

system 

of order a 


[mmy+mVniwM^io} 


[K]{x\ = p[M\[x\ 


Solve 

eigenvalue 

problem 


get pi, [ip] 


get ^ [x] 


Approximation 
to eigenvalues 


&f = pi 


* 2 

w; = Pi 


Approximation 
to eigenvectors 


ifl = mm 


[« = MM 



Tabic 2.5: Comparison of CMS and Rayleigh-Ritz Methods 
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• While the Rayleigh-Ritz .analysis uses any set of Ritz-basis vectors 
to get an approximation to the eigenvalues and eigenvectors with a 
precision that varies according to the good choice of the basis vectors, 
CMS uses, in most of the cases, the most adequate vectors (a set that 
includes the first n modes of each component), although the use of a 
truncated set precludes getting the exact solution. 

• As the number of modes used in the analysis increases, CMS converges 
monotonically from above in one of the fastest ways possible since 
each addition of a Ritz-basis vector is one of the best vectors that can 
be added. 

• In the Rayleigh-Ritz analysis, the use of a complete mass and stiffness 
matrix to get the reduced system, and then the application of a- min- 
imization process to the Rayleigh quotient is equivalent to applying 
compatibility conditions in CMS. Both methods lead to the reduced 
system of equations whose eigenvalue problem solution gives the Ritz 
coordinates [i] = p , ][V'] . 

• If the Rayleigh-Ritz analysis uses as basis vectors the eigenvectors of 
the substructures, then it leads to the same answer obtained by CMS. 

• In practical cases [K] and [M] are not known, so we do not apply 
Rayleigh-Ritz analysis, but its principle is implicitely used in CMS. 

The results obtained for the eigenvalues and eigenvectors of the assembled 
robot are approximations to the real values, and it has been shown for the 
Rayleigh-Ritz analysis [2] that the convergence is monotonic from above, 
that is, the more modes we use from the components, the closer we get 
to the real results. But this method always yields eigenvalues that are 
higher than the exact ones. The rate of convergence can be determined 
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7.00 a.oo 9.00 fO.OO 
NUMBER OF MOOES USED IN CMS. 



Figure 2.10: Convergence of 1st natural frequency with number of modes 
used 

numerically by simply comparing two sets of eigenvalues, one obtained with 
a given number of modes and the other obtained with one mode less for 
any of the components; then we can determine whether the convergence 
is satisfactory or which component needs more modes. To illustrate this 
convergence rate, consider the case of a free- free beam whose modes are 
known (see Table 2.6). If we connect it to a fixed support near the end to 
have a pinned-pinned-free beam, and we compare the first natural frequency 
of this so supported beam as obtained by CMS to that obtained from its 
formula (see Blevins [10]), we can see the effect of using more and more 
modes of the free-free beam as is shown in Figure 2.10. 

We should point out that this monotonic convergence from above can 
only be seen when the modes of the original component are correctly nor- 
malized. This can be done exactly in theory, but in practice, as we have 
discussed before (section 2.3), there is always an uncertainty on the nor- 
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malization. Therefore this convergence pattern may be different, that is, 
the convergence may be towards a set of eigenvalues other than the true 
values. 

2.7 Conditioning of the Numerical Process 

When applying CMS to connect two or more structures, and when more 
than one point of connection is used, an ill-condition may arise. Ill condi- 
tioning is a situation in which a slight change in the data produces a big 
change in the results. 

The numerical process in CMS includes the inversion of a matrix (see 
equation 2.28) plus the solution of an eigenvalue problem, and either of the 
two or both operations may be ill-conditioned. To detect this problem we 
follow the procedure described below. 

Turing [76] describes a way to determine the conditioning for matrix 
inversion; that is, given the problem [>l]{s} = {b}, we want to know 
what effect a slight change in {6} produces in the result {i}. It is found 
that Ax, = |f t A6y = (A'^ijAbj , so the biggest element in [A]~ l de- 
termines the condition number. For the case of CMS, where we have 
{p D } = [A]~ l ({Blip 1 }), we proceed as follows: 

• Run original case and get eigenvalues. 

• Find most sensitive degree of freedom in b = Bp 1 by finding j for 
^gest ^ = (A-%. 

• Alter mode shape value for most sensitive degree of freedom of B by 



1 % and get the new eigenvalues. 
Compare corresponding eigenvalues. 
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This problem of ill-conditioning was present for an early model of the robot 
joints as shown in Figure 2.11. Two of the degrees of freedom in the X- 
direction were redundant, and when connecting the robot components, the 
first resultant two modes of the robot were unexpectedly high. By applying 
the above procedure we found that the alteration of 1 % in the most sensi- 
tive degree of freedom produced 2.07 % and 1.47 % change in these two first 
modes respectively. The problem disappeared when two degrees of freedom 
in the X-direction were removed (two from the front). Then the change in 
the two modes was -0.04 % and 0.03 % and the natural frequencies were 
reasonable as expected. This elimination of two degrees of freedom in the 
X-direction did not affect the predicted robot eigenvalues appreciably, at 
least in the frequency range of interest. 

In order to identify possible redundant degrees of freedom we can pro- 
ceed as follows. In 3-D there are at least 6 degrees of freedom to be used 
since that is the minimum required to connect two rigid bodies. The way 
these degrees of freedom are distributed can not be arbitrary since they 
have to prevent any possible relative motion (translation or rotation) be- 
tween the connecting parts. If the bodies are flexible, we can start adding 
more degrees of freedom and here is where a redundant combination may 
give problems. The most sensitive degree of freedom as indicated by the 
largest element of A -1 is a good candidate at which to start looking for re- 
dundancies. In the case of the robot, an explanation of why an ill-condition 
was present is that the first two mode shapes (the ill-conditioned modes) 
are like a rigid body rocking on torsional springs at the floor, but when 
we connect the carriages on the top of the structure, we are trying to infer 
from a distant section (top of the robot) what this addition of inertia (of 
the carriages) will produce, and this addition is introduced redundantly 
since only two degrees of freedom in X-direction are necessary. 
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Y-carriage 



X-carriage Z-carriag« 



L - lead screw-nut degree of freedom 
f - friction degree of freedom 
p - rack-pinion degree of freedom 
R - pair of redundant degrees of freedom 



Figure 2.11: Redundant degrees of freedom causing ill-condition 
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Figure 2.12: Synthesized mode shape 1 of robot 

2.8 Predictions of Quasi-Static Model 

For a given configuration, that is, for a fixed position of the carriages, vi- 
bration can be described by means of the natural modes of the system. Our 
model of CMS can predict those modes and the Structural Dynamics Ana- 
lyzer can display the motions. Some of the results are shown in Figure 2.12 
and Figure 2.13. For more detail and comparison between predicted and 
measured mode shapes, see Chapter 4 and Appendix A. 
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Mode No. 


/n(Hz) 


1 





2 





3 


135 


4 


373 


5 


731 


6 


1209 


7 


1806 


8 


2522 


9 


3358 


10 


4313 



Table 2.6: Natural Modes of a Free-Free Beam 




Figure 2.13: Synthesized mode shape 2 of robot 



Chapter 3 

PROPOSED DYNAMIC 
MODEL OF A CARTESIAN 
ROBOT 



3.1 A Procedure to Analyze Systems of 
slowly-varying configuration 

Here we discuss the possibility of extending the method of CMS to systems 
which have equations with slowly varying coefficients, which is the case of 
this and many other robots as discussed below. We apply this method to 
the case of a simply supported beam with a 2-degree of freedom system 
traveling on the beam at a constant speed. Based on this simple example, 
we attempt to justify the use of this method for analysis of a cartesian robot 
and we discuss the significance of the errors resulting from the associated 
approximations. 

As discussed in Chapter 2, in the example of a moving two-degree of 
freedom system on a simply supported beam, our goal is to apply CMS 



5G 
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to connect the robot components in such a way as to take into account 
that some of the points of connection may be slowly changing with time. 
We do this by putting as the argument of the eigenfunctions of the linear 
bearing surfaces (ways) some function of time that describes how the points 
of connection are moving. This procedure, as it was applied to the beam 
problem is valid when e << 1, that is when we can neglect some terms in 
e in the general equations of motion. For a typical bridge the value of e 
falls in the range of 0.10 to 0.20. In the case of the robot, these parameters 
depend on which direction the carriages are traveling,; for each of three 
directions X,Y,Z the results are as shown in Table 3.1. Since e is small, we 
can classify the cartesian robot as a slowly-varying system; so, besides the 
fact that we can neglect some terms in e and apply CMS directly to derive 
the equations of motion, we can make further use of the slowly varying 
characteristics (if desired) to find approximate solutions to the equations 
in an efficient way. 

There is a well developed set of techniques for solving such problems 
in close form, as explained by Nayfeh [57]. The problems may be multi- 
degree of freedom and can contain forcing terms. Nevertheless, the most 
straightforward technique to solve the equations is by numerical integration, 
and that is what we have chosen to do since we first want to verify that the 
model can represent the real system with reasonable precision. 

3.2 Determination of Continuous 
Eigenfunctions 

As mentioned before, the relative motion between carriages is constrained 
by preloaded cam followers that run on linear bearing surfaces or ways, 
and physically speaking this motion is continuous. On the other hand 



CHAPTER 3. PROPOSED DYNAMIC MODEL 



58 




Figure 3.1: Linear interpolation of mode shape data 



the test points that describe the eigenfunctions of these ways are discrete 
although they are closely spaced. This situation makes it necessary either 
to interpolate between points or to curve-fit using some suitable functions. 
We first tried the interpolation techniques because they are easy to 
use. The interpolation assumes that the test points give exact values for 
the eigenfunctions; that is, that there is no error involved in the measure- 
ments. This of course is not the case. We used linear interpolation and 
cubic spline interpolation for calculating the mode shape values at the con- 
nection points, and to evaluate these techniques we looked at the resulting 
variation of the 1st natural frequency of the robot with position; this vari- 
ation should be smooth and could be checked experimentally. So, we found 
that the linear interpolation gave an irregular modal variation as shown 
in Figure 3.1, as did the cubic spline interpolation (see Figure 3.2). From 
these results we conclude that the problem must involve the presence of 
measurement errors. 
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Direction 


Dominant / n (Hz) 


Ways Length (m) 


max. speed (m/s) 


value of c 


X 


12.4 


.279 


1.27 


0.18 


Y 


18.9 


.584 


1.27 


0.06 


Z 


48.4 


.152 


1.27 


0.08 



Table 3.1: Typical Values of e for Cartesian Robot 
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Figure 3.2: Cubic spline interpolation of mode shape data 
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Next we tried curve-fitting the eigenfunctions over the test points on 
the ways, in the least squares sense, in order to smooth out the errors. 
One attempt was to use Chebyshev polynomials of different orders. These 
polynomials are defined by: 

r„(x) = 1 

T i (x)=x (3.36) 

T„n(z)=2xr B (x)-r n _ 1 
and they are orthogonal polynomials with some convenient properties for 
curve-fitting. Figure 3.3 and Figure 3.4 show the curve fits for orders 3 
and 5 for a particular case. The corresponding modal variations are shown 
in Figure 3.5 and Figure 3.6. One observation on the use of Chebyshev 
polynomials is that they give more weight to the end points of the ways 
which is an undesirable characteristic because good data may be difficult 
to get at these ends. Examination of the results showed that the modal 
variation is very sensitive to choice of the order of the polynomial fit chosen. 
This is a clear indication that the method is not robust. 

Therefore, we chose to proceed by developing a method to fit beam-like 
eigenfunctions in a least squares sense to experimental data collected on the 
ways. The ways can be considered as beams whose boundary conditions 
have to be defined. If the beam is uniform, the eigenfunction is defined by: 

4>{x) = C x sin(/3i) + C 2 cos(/?x) + C 3 sinh(/?x) + C 4 cosh(/?x) + C 5 (3.37) 

where the parameters Ci to C 5 and /? depend on the end conditions of the 
beam and x is a coordinate along the beam (see Appendix C for details). 
Curve-fitting this function to the data, in order to find the parameters, 
involves doing a nonlinear parametric fit (nonlinear in /?) which may have 
multiple solutions depending on the initial guesses for the parameters. The 
method used was again based on the Simplex technique, where we spec- 
ify any function wo want and the sum of the square of the residuals is 
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Figure 3.3: Least-squares curve-fit of mode shape data using Chebyshev 
polynomials of order 3 
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Figure 3.4: Least-squares curve-fit of mode shape data using Chebyshev 
polynomials of order 5 
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Figure 3.5: Modal variation after curve-fitting mode shape data using 
Chebyshev polynomials of order 3 
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Figure 3.6: Modal variation after curve-fitting mode shape data using 
Chebyshev polynomials of order 5 
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Figure 3.7: Two different solutions to the curve-fitting of mode shape data 
using a beam-like eigenfunction 

minimized. To find a good solution we have to do several curve-fitting 
calculations using different initial values of the parameters each time, and 
then choose one result that has the minimum standard deviation. Fig- 
ure 3.7 shows two different solutions to the curve-fitting of one case; notice 
that the best fit has a standard deviation of 5.3E-3 whereas the bad fit has 
9.6E-2. Figure 3.8 shows another case where the data is generally good and 
Figure 3.9 shows the case where there is noisy data. So, the idea of the 
equivalent beam smoothes out the errors and at the same time introduces 
the physics into the problem. When these beam-like curve fits are used to 
connect the components of the robot, the modal variation looks reasonable 
as expected (see Figure 3.10) . 

This curve-fitting technique that assumes an equivalent beam was only 
used for the cases where a motion normal to the ways would occur. For 

example, if the axis is positioned along the Y-direction, then only </>(z) and 
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Figure 3.10: Modal variation with position for case of equivalent beam 
eigenfunctions 
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4>{z) would be curve-fitted with beam-like eigenfunctions, and 4>(y) would 
be done with a straight line which would be a good model for tension or 
compression of a uniform bar. The number of curve fits done for the robot 
was as shown in Table 3.2. 

3.3 Equations of Motion 

There are four different levels of equations involved in the vibration analysis 
of the robot. These are: 

[A/]{x} + [C]{x} + [tf]{x} = {F} (3.38) 

{p} + [2^ n ]{p} + K]{p] = [4>Y{F} (3.39) 

Mtf} + WW) + [*]{«> = PTMW (3-40) 

{r} + [2 feWn ]{f} + [«2]{r} = M'PTMW (3.41) 

The first level or equation 3.38 represents the real vibrations {x} of the 
robot before it is connected; that is, when each component vibrates sepa- 
rately in an isolated way. The second level, or equation 3.39 is obtained 
from equation 3.38 by operating with the modal matrix [<f>] which decou- 
ples the equations. Equation 3.40 represents the connected robot and it 
is obtained by using the transformation matrix [T] that operates on equa- 
tion 3.39 to connect the components since it contains the compatibility 
conditions (see chapter 2). These equations are coupled again. Finally, 
level 4 is the decoupled set of equations 3.41 that is obtained from equa- 
tions 3.40 by operating with the modal matrix [if/]. To transform back to 
{x} when we know {r} we can use the following equation: 

{x} = [0][r][V]{r} (3.42) 

Before describing the procedure followed for the solution of these equations, 

we shall discuss the driving functions. 
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3.4 Force Analysis and Prescribed Motion 

So far in this thesis we have dealt with the left-hand side of the equations, 
that is, the part that describes the homogeneous solution of the system. The 
right-hand side which contains the forcing terms that excite the system 
are described now. When the carriages move, the inertia forces and the 
gravity forces acting on a system that changes configuration tend to excite 
vibrations on the structure, and so the forcing terms are associated with 
the motion of the carriages. 

In order to leave out the control system that causes the motion, and 
consider only the structural vibrations, we measure the motion of the motor 
rotor including torsional vibration during a prescribed robot motion and 
use this information as prescribed motion at the motor rotor to drive the 
equations of motion. 

With regard to the interface between the carriages and structure in the 
model we can do the simulations in two different ways: 

1. At a given instant the carriages can be assumed not have relative 
motion, so they are pinned to the ways, but those pins are slowly 
changing with time. Here we subtract the rigid body mode by as- 
suming that the carriages are rigid and calculating the reaction forces 
that occur due to an acceleration equal to the rigid body acceleration. 
These reaction forces are then used to drive the equations of motion. 

2. The carriages can be assumed not to be pinned to the ways so that 
they can slide in the direction of motion. However they have to over- 
come a friction force which acts as an external force whose value we 
know. Then, at the end of the move we set this force to zero. 

Since this second alternative does not require knowledge of the rigid 
body acceleiation, which strictly speaking is not known, wc have chosen to 
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do the simulations of the robot vibrations using this second method. 

3.4.1 Prescribed Motor-rotor Vibration 

In this procedure we eliminate the forces exerted by the motors from the 
equations of motion in favor of the prescribed vibration which acts as a 
"forcing" function to drive the system. 

The problem, besides being position-dependent presents limitations in 
the available data, that is, all we know are the natural modes of each 
component and a series of complex frequency responses from which the 
modes were calculated. The problem may be solved in several ways as 
explained in Appendix B, but the one that we consider the most suitable 
for the characteristics of this cartesian robot is a time-domain technique 
which is not restricted to time- invariant systems. 

As we have discussed in Chapter 2, the application of the compatibility 
conditions in the CMS method leads to a set of coupled differential equa- 
tions (equation 2.11 ). These equations represent the complete cartesian 
robot, and they are driven by forces exerted by the motors in the presence 
of friction and gravity forces. These equations can be decoupled first for 
convenience and then manipulated to eliminate the motor forces by sim- 
ple multiplications by constants and adding together different rows of the 
system, since the motor forces appear in linear form. In doing this elimi- 
nation, the equations become coupled again, and we reduce the number of 
equations by the number of motor forces eliminated (maximum of 3 forces, 
one for MX, one for MY, and one for MZ). Consequently we must add new 
equations, and these come from the relations between the variables r of 
the assembled system and the real encoder positions Imx> x my, or x mz as 
follows 

XMX = l<t>Mx}[T][1>]{r},etc. (3.43) 
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where [^x] is a row vector coming from the modal matrix evaluated at 
the motor shaft. 

So, by adding these equations we again have the number of equations 
required for the number of unknowns ({r}). The resulting equations can 
be represented by 
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(3.44) 
where f c contains the system external forces other than motor forces. These 
equations can now be integrated directly to find {r}, and from there the 
vibration at any point i can be obtained from {x} = [0][T]ty]{r}. 

Notice that we have differentiated equation 3.43 twice with respect to 
time in order to obtain a mass matrix which remains non-singular. This is 
a convenient situation to have for numerical integration of the equations. 
This differentiation of equation 3.43 assumes that the terms in 4> and <f> are 
negligible compared to the terms in f and f respectively which is correct 
for a slowly-varying system. For the vibration analysis of the moving robot 
we moved it always in one single direction (X,Y, or Z) to simplify the 
mathematical manipulations. 



3.4.2 Friction Considerations 

The friction force plays an important role in the robot vibrations. This 
force is high due to preload in the cam followers; it normally represents 
about 15 percent of the total motor force during motion. 

The friction force has been measured by Benjamin [5] who did a computer- 
controlled experiment on the cartesian robot to determine the values of the 
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Figure 3.11: Variation of friction force along x-direction 



friction force along the X-direction of the robot for different positions of 
the X-carriage. He modeled this force variation by a straight line which 
was fitted to a number of points obtained experimentally. His results are 
reproduced in Figure 3.11. For the Y-direction and Z-direction he did not 
report any results. 

A simple static test measuring the force required to start motion was 
done to get approximate values of the friction force for four different posi- 
tions along the ways for each direction. From these values we approximate 
a position-dependent function that describes the friction force for any of 
the three principal directions. 

These forces are used in the model in an attempt to simulate what 
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physically happens during the motion. We assume that the friction force is 
equal to the saturation value of friction (follows the straight line variation) 
as long as there is relative motion between carriages and ways. It may 
happen that the friction forces reverse sign, and that is when the difference 
in displacement between two sliding parts changes sign. This occurs for the 
case of motion in x-direction when : [<f> x ~ carriage] [T] [V>] - [^y- carriage] [7"] [V»] 
changes sign, where <t>x -carriage is the mode shape value at the instantaneous 
slide point of the X-carriage, and (f>y - carriage is the corresponding value of 
the Y-carriage. The program has to monitor this difference and adjust the 
sign of the friction force accordingly. 

As soon as the motion stops, the friction force is assumed to be always 
fluctuating between the— saturation value and the + saturation value while 
not allowing any relative motion. Under these circumstances the force is 
no longer known. To cope with this we pin the carriage to the ways using 
friction pins. This means that we modify the model by connecting the 
carriages to the ways at the end of a simulated move. This gives a slightly 
different dynamic behavior of the structure. In practice we noticed better 
correlation with the experimental results when we simply set the friction 
force to a constant value during the move and zero at the end of the move; 
this suggests that perhaps there is a little backlash in the nut-lead screw 
connection. 

3.4.3 Conditioning of the Encoder Motion 

The encoder motion is used as a prescribed motion at the motor rotor 
to drive the equations of motion. The optical encoder used gives a dis- 
placement (rotation) function which is stepwise continuous as shown in 
Figure 3.12. This displacement when used without any conditioning to 
calculate velocity and acceleration gives strong discontinuities mainly in 
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Figure 3.12: Encoder displacement 

the acceleration function because the encoder provides discrete rather than 
continuous information (see Figure 3.13 and Figure 3.14). To avoid these 
problems we have to filter out the stepwise characteristic of the encoder 
measurements; one way is by using a digital filter. 

The filter considered is a differentiating, low-pass filter with zero phase- 
shift which besides differentiating the encoder displacement, filters out the 
undesired higher harmonics. 

This digital filter is a non-recursive type filter because the filter output 
g m is calculated explicitly from the input f m by means of 



N 



9m — / > "njm-n 

n=-AT 



(3.45) 



where b n are the filter coefficients. This summation is an approximation to 
the convolution integral where the set of 6 n constitutes the impulse response 
of the filter. These coefficients have to be calculated so that the frequency- 
domain effects of this process come out as desired, that is, the desired 
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Component 


No. of ways 


No. of modes 


No. of fits 


STR 


2 


13 


78 


Y 


2 


10 


60 


Z 


4 


5 


60 



Table 3.2: Number of Curve Fits for Ways 




Figure 3.13: Unfiltcred encoder velocity 
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Figure 3.14: Unfiltered encoder acceleration 



b n =£- [ T iy d (y w )e jnwr dw, -N<n<N 



(3.46) 



transfer function of the filter Hd(ju) can be approximated by some H[jui) 
through a least squares process. In this way we find 6 n for n in the interval 
i-N,N}. 

The general formula for the coefficients, as given by Stearns [72] is 

T 
l-r 

We chose to design a filter which would eliminate all frequencies above 
130 Hz (this limit being set by the first 10 modes of the robot), and differ- 
entiate the displacement history twice as well. An algorithm developed by 
McClellan, Parks and Rabiner [51] can do this linear phase filter design, al- 
though its differentiation is of first order, so we have to do the process twice 
starting with displacement to get acceleration. The number of points(N) 
was 63. 

The impulse response obtained for this filter is shown in Figure 3.15 
and the corresponding Fourier transform is shown in Figure 3.16, whore 
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Figure 3.15: Impulse response of digital filter 

the Nyquist frequency is F = 1500 Hz which corresponds to the sampling 
time t = 0.67 ms of the computer-controlled experiment. 

The application of the digital filter to the encoder data gives the velocity 
and acceleration shown in Figure 3.17 and Figure 3.18 respectively. These 
data were collected during a move in the X-direction. Note the component 
of acceleration at 87.5 Hz which corresponds to the dominant 8th mode of 
vibration (see Appendix A). 



3.5 Solution of the Equations of Motion 

In order to calculate the time history x of vibrations at any of the test points 
of the robot we can integrate equation 3.41, truncating to a convenient 
number of modes if we want, and then using equation 3.42 we can transform 
back to real coordinates. This transformation has to be done at each step 
since the matrices involved in the calculation are position-dependent. Now, 
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Figure 3.16: Fourier transform of digital filter 
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Figure 3.17: Filtered encoder velocity 
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Figure 3.18: Filtered encoder acceleration 

if we want to have a fast computer model to integrate these equations, we 
have to do all the work that can be done prior to the integration, so that 
the integration itself is fast and can hopefully be done in real time. For this 
purpose, three computer programs were prepared to solve the equations as 
described below: 

• The first program calculates the modal parameters at discrete points 
along the path, say at every 10 mm. 

• The second program smoothes the modal parameters using cubic 
splines to calculate these parameters at each time step. These two 
programs can be run prior to the integration. 

• The last program integrates the equations by using the parameters 
derived above at each time step. It integrates equation 3.41 with a 
4th order Runge-Kutta-Merson algorithm. Then it transforms r to x 
to generate a time history at the desirod test points. 



CHAPTER 3. PROPOSED DYNAMIC MODEL 79 

This scheme was implemented on a VAX 11/750 computer and the ratio 
of computation time to real time was 40 to 1 when we integrated over the 
first 10 modes of equation 3.41. To improve this speed for real time control 
purposes we could proceed in a number of ways, for example using a faster 
computer or trying to use fewer equations in equation 3.41, or even trying 
to get asymptotic approximations to these equations with slowly-varying 
coefficients, but this is beyond the scope of the present work. 

3.6 Some Results of the Dynamic Model 

To test the model, a series of experiments were conducted on the robot 
which are explained in detail in chapter 4. From these tests, the motion 
at the motor-rotor for the cases of motion along the X-direction, or Y- 
direction, or Z-direction could be determined, and it is from this data that 
we can prescribe the motion of the carriages to run the dynamic model. 
Figure 3.19 shows the vibration at the end-point (test point 144) of the 
robot when the carriages are moved in the X-direction. This curve is a time 
domain representation and covers the whole motion, that is from t = to 
t = 0.57 s. Figure 3.20 shows the robot tip residual vibration (t = 0.57 to 
t = 0.89 s). Figure 3.21 shows the frequency content of this end vibration. 
For the case of the Y-direction and Z-direction see Chapter 4 where we 
compare these simulation results to the test results. 
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Figure 3.19: Vibration of end-effector in x-direction 




Figure 3.20: End vibration in x-direction 
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Figure 3.21: Fourier transform of end vibration in x-direction 



Chapter 4 

EXPERIMENTS AND 
RESULTS 



4.1 Objectives 

We are trying to understand how a time-dependent system vibrates in order 
to be able to control the vibration and/or design the system in a better way. 
For this purpose we have developed a mathematical model which must then 
be validated by test. 

There are two classes of objectives for the tests described in this thesis, 
one is that we need modal data from each of the robot components to be 
used as input to the mathematical model. The other class of objectives 
is the mathematical model validation which should evaluate the accuracy 
of this model. The modal tests which were done to gather input data 
are briefly mentioned in this chapter, but fully described in Appendix A. 
In contrast, the validation tests are explained with more detail here, and 
a comparison with the predicted results from the mathematical model is 
presented. 

Wc will discuss four different types of tests. The first two provided data 
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for the model; the last two were used in verification. The tests included 

• Modal tests to determine natural frequencies, modal damping and 
mode shapes. 

• High- resolution transfer function measurements for orthonormaliza- 
tion of eigenvectors. 

• Measurements of modal variation with position of the carriages. 

• Measurements of robot vibration when the carriages are moving. 

4.2 Modal Tests 

The theory of modal analysis deals with two kinds of modes, normal modes 
and complex modes. The type of modes depends on the damping behav- 
ior of the structure. Caughey [14] presents an analysis of the conditions 
under which a damped linear system possesses classical normal modes. He 
concludes that the damping matrix must be diagonalized by the same trans- 
formation that uncouples the undamped systems. Meirovitch [53] gives a 
modal analysis procedure designed to handle the general case of damping. 
The modal analysis as is done in the laboratory with the use of digital 
equipment may utilize either of the two methods, although in general real 
structures should be analyzed using complex mode theory. Ibrahim [40] 
explains the errors involved in the normal mode approximation to complex 
modes, concluding that even for lightly damped structures this approxima- 
tion may lead to large errors. Potter [62] and Potter and Richardson [63] 
explain the theory of complex modal analysis showing how to calculate 
modal vectors from a measured transfer matrix. 

With regard to vibration tests, Mustain [56] shows the results of a survey 
of modal vibration test/analysis techniques in the USA. Favour, Mitchell, 
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and Olson [23] discuss the historical development of transient test tech- 
niques utilizing a digital computer. Halvorsen and Brown [28] describe the 
impact procedure, its theory, applications and limitations. Ramsey and 
Richardson [65] present a brief review of modal analysis based upon the 
identification of modal parameters from measured transfer functions. They 
point out the advantages of a band selectable Fourier transform (zoom 
transform) for obtaining increased accuracy, resolution and dynamic range 
in the measurements. Richardson and Kniskern [66] present an algorithm 
that reduces error in the mode vector when using more than one simple row 
or column of residue data( more than one excitation point when using a 
shaker, or more than one response point when using a hammer). Okubo [59] 
analyzes the effect of nonlinearities on transfer function measurements with 
impact excitations. With regard to nonlinearities in Structural Dynamics, 
Crandall [19] discusses the nature of nonlinear structural models and sur- 
veys a variety of nonlinear response phenomena which can be predicted by 
such models. 

For the modal tests, we have in the Acoustics and Vibration Laboratory 
an HP-5423A Structural Dynamics Analyzer. In conjunction with measur- 
ing instruments to determine acceleration and force, and with excitation 
means like an electrodynamic shaker or an instrumented hammer, the An- 
alyzer provides adequate facilities for our purposes. This Analyzer assumes 
proportional damping when doing modal analysis calculations. 

The vibration tests performed on the robot components to obtain input 
data for the mathematical model are explained in detail in Appendix A. 
Basically we obtain from 4 to 10 elastic modes for each component, which in 
addition to the corresponding 6 rigid body modes, give a sufficient number 
of modes for the model to provide convergence. 

Of more interest in this chapter is the validation test performed on the 
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completely assembled robot which was tested in the configuration shown in 
Figure 4.1. Here the X-carriage and Y-carriage are in the middle of their 
respective way paths, and the Z-carriage is all the way down. The test 
set-up is shown in Figure 4.2. A Bruel and Kjaer electrodynamic shaker 
type 4801 with mode study head type 4814 was used to excite the structure 
at the test point no. 145 in the -Z direction while a triaxial accelerometer 
B & K type 4340 was used to transmit acceleration to the analyzer via 
a charge preamplifier B & K type 2628. The force was measured by a 
Wilcoxon impedance head attached between the shaker and the structure. 
For each of the test points where the accelerometer was attached, and for 
each direction (X,Y, or Z) a transfer function made of 50 to 80 averages 
was calculated and stored in the analyzer. The frequency range was 0-200 
Hz in the first set of measurements, and 0-25 Hz in the second set with a 
total number of measurements of 405 for each test. 

A typical transfer function is shown in Figure 4.3, and the natural fre- 
quencies and modal damping factors are shown in table 4.1. Here in this 
table we also show the corresponding predicted values obtained by the 
mathematical model. The first 4 mode shapes as obtained by test are 
shown in the top part of Figures 4.4- 4.7 whereas the corresponding pre- 
dicted modes are shown in the bottom part. 

Comparing the experimental and synthesized results we notice that the 
average error in the first 17 natural frequencies is 6.9 %. The mode shapes 
show the same type of tridimensional motion for both sets of results. 

We believe that this agreement is close enough and that the math- 
ematical model can also predict the modal characteristics for any other 
configuration with comparable accuracy. 
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Figure 4.2: Test set-up for modal analysis 




Figure 4.3: Typical transfer function 
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Mode 
No. 


Meas 

/» (HZ) 


ured 

in (%) 


Predicted 
/ n (Hz) in (%) 


1 


12.3G 


1.98 


11.93 


1.45 


2 


18.85 


1.65 


19.17 


1.45 


3 


36.36 


1.42 


36.87 


0.78 


4 


42.83 


9.52 


46.34 


3.30 


5 


47.83 


4.90 


49.53 


2.54 


6 


60.51 


3.24 


61.38 


1.45 


7 


70.72 


0.89 


75.06 


1.82 


8 


81.47 


1.95 


96.47 


0.63 


9 


86.67 


1.74 


101.75 


1.41 


10 


97.70 


1.17 


123.39 


1.30 


11 


112.52 


7.39 


126.84 


1.06 


12 


126.27 


0.10 


128.72 


0.75 


13 


139.43 


2.62 


141.05 


1.16 


14 


149.85 


1.00 


145.68 


2.86 


15 


152.34 


0.66 


148.23 


2.76 


16 


178.12 


1.79 


169.39 


1.04 


17 


184.37 

■ ■ —J , 


0.01 


179.06 


1.11 



Table 4.1: Natural Frequencies and Modal Damping Factors. 
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Figure 4.4: Mode shape no. I of robot as obtained by test (top) and by 
synthesis (bottom) 
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Figure 4.5: Mode shape no. 2 of robot as obtained by test (top) and by 
synthesis (bottom) 
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Figure 4.6: Mode shape no.3 of robot as obtained by test (top) and by 
synthesis (bottom) 
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Figure 4.7: Mode shape no. 4 of robot as obtained by test (top) and by 
synthesis (bottom) 
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4.3 High-Resolution Transfer Functions 

The HP analyzer calculates transfer functions using 256 complex numbers 
for equally spaced frequencies in the specified frequency bandwidth. This 
resolution can be increased by repeating the test with different center- 
frequencies, and so, ending up with several transfer functions that can later 
be collected into one single plot outside of the analyzer, in a computer. 
We transmitted these transfer functions to a PDP11/44 computer through 
a GPIB interface. The resultant high-resolution transfer function which 
contains detailed information around each mode of interest can be used 
to get more accurate estimates of the modal parameters, particularly the 
mode shape value for orthonormal eigenvector as explained in chapter 2. 

Two of the robot components were tested for these high-resolution trans- 
fer functions; the robot structure (without the carriages) attached to the 
floor, and the Z-carriage. Details of the test are given in Appendix A. 
One of these curve-fits for the Z-carriage is shown in Figure 4.8. As was 
discussed in chapter 2, the results may seem accurate but there is an impor- 
tant variance when we apply different curve-fitting methods. For example, 
if we change the initial conditions of this nonlinear parametric fit we get 
the curve fit shown in Figure 4.9 which may seem a very similar result but 
the mode shape values differ from 5 to 25 percent depending on the mode, 
and that happens because damping changes as well in such a way that the 
combination of <f>kk, <^k, and ft gives another solution of the nonlinear prob- 
lem. Therefore, this difficulty in the system parameter identification may 
preclude good accuracy of the mathematical model. 
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Figure 4.8: Multimodal curve fit no.l of transfer function of Z-carriage 
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Figure 4.9: Multimodal curve fit no.2 of transfer function of Z-carriage 
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4.4 Modal Variation with Position 

A position-dependent system changes its modal characteristics as it changes 
configuration. The cartesian robot shows slight changes when the rela- 
tive position of the carriages change. In order to partially evaluate these 
changes, the robot was tested in different configurations obtaining for each 
one a single transfer function which may show at least the variation in 
natural frequencies. Table 4.2 shows these variation for the 1st mode as 
well as the corresponding calculated value from the mathematical model. 
Table 4.3 shows similar results for the second mode. 

These tables (4.2 and 4.3) and the set of mode shape plots show in gen- 
eral good agreement between experimental and calculated results. Based 
on this and on the results of section 4.1 we can generalize the validity of 
the model to the whole work space of the cartesian robot. Table 4.4 shows 
the variation of the first ten natural frequencies along the Y-direction for 
different positions in the X and Z directions (calculated from the model). 
The values of Y (mm) represent the position of the Y-carriage as measured 
in the global system of coordinates. 

4.5 Vibration of Moving Robot 

The final objective of this thesis is to be able to predict the vibrations of 
any part of the robot as the carriages move, and the experimental results 
described in this section permit evaluation of the model simulations. Three 
different experiments were carried out: 

• Motion of the carriages in the X-direction only. 

• Motion in Y-direction. 

• Motion in Z-dircction. 



CHAPTER 4. EXPERIMENTS AND RESULTS 



96 



Position of carriages 


fi (Hz) measured 


/i (Hz) predicted) 


■A-fronti * middle i '"down 


11.80 


11.98 


■^ front > J/p/ti ^(ioton 


11.22 


11.70 


■"■middle* * middle i '"down 


12.36 


12.05 


^middle i 'left i Zdown 


10.94 


11.75 


■"■midtilm * middle i "up 


11.49 


12.27 



Table 4.2: Modal Variation with Position, first Mode. 



Position of carriages 


/2 (Hz) measured 


fi (Hz) predicted) 


■"■ front i * middle ! '"down 


18.90 


19.09 


X front, Ylefti "down 


19.03 


19.23 


■X-middlei * middle > "down 


18.85 


19.17 


X-middle-, Hcftt "down 


17.97 


19.31 


■''■middle > * middle i ^up 


18.82 


19.15 



Table 4.3: Modal Variation with Position, second Mode. 
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Table 4.4: Variation of Natural Frequencies (Hz) with Position for the First 

Ten Modes, according to the model. 
X back Z top 

Y= 422.28: 12.1 19.3 3G.9 47.9 49.9 61.1 69.9 88.7 101.9 105.4 
Y= 578.07: 12.4 19.2 36.9 48.3 49.8 61.6 71.0 91.0 102.1 110.4 
Y= 733.85: 12.5 19.1 36.9 47.9 49.9 61.8 72.0 92.5 101.1 106.8 
Y= 889.64: 12.7 19.1 36.9 47.6 50.0 62.0 72.5 93.1 100.4 107.8 
Y=1006.48: 12.9 19.0 36.9 47.2 50.1 62.1 72.9 93.5 99.8 108.1 
X bark Z middle 

Y= 422.28: 12.1 19.3 36.9 47.7 49.9 61.2 69.5 90.1 101.6 119.2 
Y= 578.07: 12.4 19.2 36.9 48.0 49.9 61.7 70.8 92.7 106.1 121.5 
Y= 733.85: 12.6 19.1 36.9 47.6 50.0 61.9 71.8 93.7 101.8 122.9 
Y= 889.64: 12.9 19.1 36.9 47.3 50.2 62.1 72.7 94.2 102.0 123.4 
Y=1006.48: 13.0 19.0 36.9 46.9 50.3 62.2 73.2 94.3 101.7 122.5 
X back Z down 

Y= 422.28: 11.6 19.3 36.9 47.5 49.4 61.0 67.3 89.5 99.7 121.3 
Y= 578.07: 11.8 19.2 36.9 47.8 49.6 61.4 68.7 91.8 103.9 124.2 
Y= 733.85: 11.9 19.1 36.9 47.4 49.7 61.6 69.7 92.2 100.5 126.2 
Y= 889.64: 12.1 19.1 36.9 47.0 49.9 61.7 70.5 92.4 100.7 126.9 
Y=1006.48: 12.2 19.0 36.9 46.6 50.0 61.8 71.1 92.3 100.4 125.6 
X middle Z top 

Y= 422.28: 11.9 19.3 36.9 47.3 49.2 60.9 74.7 89.5 100.6 118.0 
Y= 578.07: 12.2 19.2 36.9 47.4 49.3 61.4 75.9 92.3 101.2 121.1 
Y= 733.85: 12.3 19.1 36.9 47.0 49.4 61.6 76.5 94.3 101.6 120.5 
Y= 889.64: 12.4 19.1 36.9 46.6 49.5 61.7 76.9 95.4 101.4 119.6 
Y=1006.48: 12.5 19.0 36.9 46.2 49.6 61.8 77.1 95.6 101.0 119.8 
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Table 4.4: continued 
X middle Z middle 

Y= 422.28: 12.2 19.3 3G.9 47.2 49.5 61.1 74.1 92.2 100.7 121.8 

Y= 578.07: 12.5 19.2 30.9 47.3 49.6 61.5 75.5 95.4 102.6 122.3 

Y= 733.85: 12.6 19.2 36.9 47.0 49.7 61.7 76.4 97.0 102.3 122.3 

Y= 889.64: 12.9 19.1 36.9 46.5 49.8 61.9 77.0 97.5 102.7 121.9 

Y=1006.48: 13.0 19.0 36.9 46.1 49.9 62.0 77.4 97.1 103.1 121.6 

X middle Z down 

Y= 422.28: 11.8 19.3 36.9 46.9 49.2 60.9 73.0 92.4 100.9 121.3 

Y= 578.07: 11.9 19.2 36.9 46.9 49.3 61.3 74.4 95.0 103.3 122.3 

Y= 733.85: 12.0 19.2 36.9 46.5 49.5 61.5 75.4 96.4 103.0 122.0 

Y= 889.64: 12.2 19.1 36.9 46.1 49.6 61.6 76.1 96.1 102.9 121.6 

Y=1006.48: 12.3 19.0 36.9 45.7 49.6 61.6 76.5 95.4 103.2 121.4 

X front Z top 

Y= 422.28: 11.5 19.2 36.8 46.1 48.7 60.5 75.0 87.6 101.2 117.2 

Y- 578.07: 11.6 19.1 36.9 46.1 49.0 60.9 76.1 90.8 101.9 120.5 

Y= 733.85: 11.7 19.1 36.9 45.8 49.1 61.1 76.3 93.1 102.4 121.9 

Y-- 889.64: 11.7 19.0 36.8 45.3 49.2 61.1 76.7 93.9 102.1 123.1 

Y=1006.48: 11.8 18.9 36.8 44.9 49.2 61.2 76.9 94.0 101.7 123.4 

X front Z middle 

Y= 422.28: 11.9 19.2 36.9 46.3 48.9 60.8 75.4 89.9 102.2 115.2 

Y= 578.07: 12.1 19.1 36.9 46.2 49.1 61.2 76.8 93.8 102.8 119.1 

Y= 733.85: 12.2 19.1 36.9 45.9 49.2 61.3 77.0 96.2 102.6 120.1 

Y= 889.64: 12.4 19.0 36.9 45.5 49.3 61.4 77.4 97.2 102.8 121.3 

Y=1006.48: 12.5 19.0 36.8 45.1 49.4 61.5 77.6 96.9 102.9 121.7 

X front Z down 

Y= 422.28: 11.7 19.2 36.8 46.1 48.7 60.6 75.1 89.4 92.1 113.3 

Y= 578.07: 11.9 19.2 36.8 46.0 48.9 61.1 76.6 89.4 95.7 116.7 

Y= 733.85: 12.0 19.1 30.8 45.7 49.1 61.2 76.8 89.2 97.8 117.4 

Y= 889.64: 12.1 19.0 30.8 45.2 49.2 01.3 77.4 88.2 99.2 118.3 

Y HU0G.I8: 12.2 l<>.0 36.8 1 1.8 19.2 61..'! 77. 6 87.5 99.0 118.5 
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Each test gives out two results, the vibration of the end-effector as mea- 
sured with an accelerometer, and the motor-rotor displacement as measured 
with an optical encoder. 

The dynamic model which was discussed in the previous chapter may 
now use this motor motion as input to simulate the robot vibrations, and 
therefore an evaluation can be done by comparing the experimental to the 
simulated results. This comparison is shown in detail below. 

4.5.1 Computer- Controlled Tests 

In order to get information about the robot vibration as the carriages are 
moving, we made use of the PDP11/23 computer which controls the robot 
to interact with the rest of the instrumentation as depicted in Figure 4.10. 
This set-up is the same used by Singer [71]. The motion specified to the 
carriage motor (one at a time) is controlled by a PWM pulse width modu- 
lated amplifier which is driven by the computer trying to follow a velocity 
profile (linearly increasing during the first half and linearly decreasing in 
the second half). The amplifier drives the motor and receives feedback from 
it to correct for deviations to the specified motion. The real velocity profile 
( Figure 4.11) for a test along X-direction shows the corrections done when 
compared with the specified profile. 

The torsional rigid body motion and superimposed vibration for any of 
the three motors MX, MY, or MZ are measured by an incremental optical 
encoder HP type HEDS-6010 B08 which transmits data to a microcom- 
puter connected to the main computer. These data can then be transmitted 
to a bigger computer (VAXl 1/750) to form a data file. 

The acceleration at the end-point (or any other of the test points) is 
received by the HP analyzer in the usual way, but the trigger is again 
controlled by the computer through a D/A converter. 
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Figure 4.10: Test set-up for computer-controlled experiment 
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Figure 4.11: Actual Velocity Profile followed by Controller. 

4.5.2 Vibration Results for Motion in X-direction 



With the Y-carriage in the middle and the Z-carriage down we ran the MX 
motor along the X-direction from Xi = 154 mm to X 2 = 291 mm with an 
average acceleration of 0.176ff and an average deceleration of O.I8I3. The 
resulting motor displacement is shown in Figure 4.12, and the corresponding 
velocity after being filtered is given in Figure 4.13. 

The end-point vibration as measured by an accelerometer pointing along 
the X-direction is shown in the top part of Figures 4.14 and 4.15 (complete 
vibration and vibration at the end of the move, respectively). The linear 
spectrum of the end vibration is shown in the top part of Figure 4.16. The 
corresponding results as obtained with the mathematical model are shown 
in the bottom part of the same figures. 
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Figure 4.12: Encoder displacement during motion in X-direction 
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Figure 4.13: Encoder velocity during motion in X-direction 
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Figure 4.14: Vibration of end-effector during motion in X-direction as ob- 
tained by test (top) and by simulation (bottom) 
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Figure 4.15: Vibration of end-effector at the end of the move for motion in 
X-dircction as obtained by test (top) and by simulation (bottom) 
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Figure 4. 16: Fourier transform of vibration of end-effector at the end of the 
move for motion in X-direction as obtained by test (top) and by simulation 

(bottom) 
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Figure 4.17: Encoder displacement during motion in Y-direction 
4.5.3 Vibration Results for Motion in Y-direction 

The X-carriage was positioned in the middle and the Z-carriage down when 
the carriages were run in the Y-direction from Vi = 498 mm to K 2 = 683 
mm with an average acceleration of 0.25g, and an average deceleration 
of 0.32g. The motor displacement and velocity are shown in Figure 4.17 
and Figure 4.18 respectively. The vibration results as obtained by test are 
shown in the top part of Figures 4.19 to 4.21. The corresponding results 
for the model are shown in the bottom part. 



4.5.4 Vibration Results for Motion in Z-direction 

With the X-carriage and Y-carriage positioned in the middle we ran the 
Z-carriage from Z\ — -486 mm to Z-i — -372 mm with an average accelera- 
tion of 1.51(7, and an average deceleration of 0.76g. The motor displacement 
and velocity are shown in Figure 4.22 and Figure 4.23 respectively. Here 
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Figure 4.18: Encoder velocity during motion in Y-direction 

the encoder displacement was filtered using a cut-off frequency of 270 Hz 
because the test was of shorter duration and the analyzer could sample the 
resulting vibration with a larger frequency bandwidth. The vibration re- 
sults as obtained by test are shown in the top part of Figures 4.24 to 4.26. 
The corresponding results of the model simulation are shown in the bottom 
part. 



4.5.5 Comparison of Simulation and Experimental Re- 
sults 

The final test of the mathematical model is in these dynamic experiments 
where the robot carriages are run in all three directions and several mea- 
surements are taken at the robot end-effector. A careful study of these 
results as compared to the model results gives the following general obser- 
vations. 
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Figure 4.19: Vibration of end-effector during motion in Y-direction as ob- 
tained by test (top) and by simulation (bottom) 
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Figure 4.20: Vibration of end-effector at the end of the move for motion in 
Y-direction as obtained by test (top) and by simulation (bottom) 
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Figure 4.21: Fourier transform of vibration of end-effector at the end of the 
move for motion in Y-direction as obtained by test (top) and by simulation 

(bottom) 
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Figure 4.22: Encoder displacement during motion in Z-direction 




Figure 4.23: Encoder velocity during motion in Z-dircction 
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Figure 4.24: Vibration of end-effector during motion in Z-direction as ob- 
tained by test (top) and by simulation (bottom) 
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Figure 4.25: Vibration of end-effector at the end of the move for motion in 
Z-direction as obtained by test (top) and by simulation (bottom) 
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Figure 4.26: Fourier transform of vibration of end-effector at the end of the 
move for motion in Z-direction as obtained by test (top) and by simulation 

(bottom) 
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• The model predicts the same order of magnitude of the results and 
the same dominant harmonics in the spectra. 

• The encoder acceleration, after being filtered, contains the dominant 
harmonics which are excited by the motor, harmonics which in turn 
reflect back to the excitation when the motor tries to maintain a given 
velocity profile. 

• The system responds not only with the first two or three modes, but 
rather with a broad band of frequencies which may be excited. 

• The vibration that occurs during the move, in the direction of motion, 
follows very closely that of the motor rotor. The end vibration how- 
ever, has its own characteristics which are in part due to the motor 
whose forces remain active in order to fix the rotor. 

A more specific discussion about each of the motion directions is given now: 

Motion in X-direction 

• The vibration during the move and the end vibration as measured 
by the accelerometer shows some drift which if corrected would make 
the plot of experimental data look more like the calculated vibration. 

• The end vibration has a peak-to-peak amplitude of approximately 
0.30 g for the test case and 0.24 g for the calculated one. 

• There are two dominant frequencies, one at 12.3 Hz and another at 
88 Hz which correspond to the robot modes no.l and no. 8 as shown 
in Figures ?? and ??. These modes are clearly excited because of 
their predominant X-direction component of end point motion. 
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• The relative amplitude of these harmonics as seen in the frequency 
spectra of Figure 4.16 tend to agree for both cases, experimental and 
simulated. 

Motion in Y-direction 

• The vibration during the move shows the same general shape for 
both cases with a peak-to-peak envelope amplitude of 2.0 g for the 
test case, and 2.5 g for the model. 

• The end vibration shows for the test case an envelope of 0.85 g, and 
for the model case an envelope of 1.2 g. 

• The dominant frequencies for both cases have the same general ten- 
dency as shown in Figures 4.19 and 4.20. The two most important 
harmonics occur at 19 Hz and 36 Hz, and correspond to the robot 
modes no. 2 and no.3. 

Motion in Z-direction 

• The vibration during the move has an envelope of 5.0 g for the test 
case, and 7.0 g for the model case, the difference being that the model 
responds more at the higher dominant harmonic of 230 Hz. 

• The end vibration as obtained by test shows more drift here than 
in the previous tests. Damping decay is also more visible in this 
direction. 

• The first two dominant frequencies in the end vibration are 115 Hz 
and 176 Hz (perhaps multiples of the line frequency, 60 Hz), but 
in general all harmonics are rather uniform in amplitude as seen in 
Figure 4.26. The spectra in both cases look similar. 



Chapter 5 

CONCLUSIONS AND 
RECOMMENDATIONS 



5.1 Conclusions 

1. The basic idea behind the present vibration analysis of the cartesian 
robot is the connection of substructures which remain time-invariant, 
and whose modal characteristics and transfer functions can be read- 
ily determined either experimentally or analytically. For analytically 
determined data, we normally use a finite element analysis proce- 
dure to connect the substructures, but CMS can also be used. For 
experimentally determined data, the classical way of connecting sub- 
structures is the Impedance method, although CMS can be used here 
also. For this particular problem, the position dependency eliminates 
the possibility of application of the Impedance method, and for the 
slowly-varying case CMS and FE can be used. The relatively com- 
plex structure of each component and the physical availability of the 
prototype point towards the use of experimental data for getting a 
more realistic model. Therefore, the use of CMS becomes logical. 
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2. The fact that CMS uses modal data as input permits us to process 
all experimental data separately, and input the processed data to the 
model which can then perform fewer operations to predict the robot 
vibrations in a fast way. The matrix handling features of CMS are 
very convenient when a relatively large number of degrees of freedom 
are to be used to connect the substructures, as in the case of the 
robot. 

3. The extension of CMS to the case of slowly-varying systems can be 
done when the speed at which the carriages move is such that the 
time spent in traveling all along the ways is much bigger than a typ- 
ical natural period of oscillation. This situation allows us to neglect 
terms in <j> and <j> and permits a simple manipulation of the data. 
The example of a traveling 2- degree of freedom system on a simply 
supported beam analyzed with two time scales in Chapter 3 gives an 
illustration of the validity of these approximations. 

4. When connecting substructures through multiple degrees of freedom 
it may happen that a singularity or an ill-condition is present. Since 
the CMS process involves both the inversion of a matrix A that incor- 
porates the compatibility conditions, and the solution of an eigenvalue 
problem, a bad distribution of degrees of freedom causes the eigenval- 
ues to vary significantly from expected values. Furthermore a slight 
change in the data may cause a large variation in some eigenvalues. 
This can be avoided by detecting those redundant degrees of freedom 
and eliminating them or changing the distribution. A technique has 
been developed in this Thesis to detect and solve these problems. 

5. Modal parameter identification has been an issue when we want high 
precision evaluation of the modal parameters based on experimen- 
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tal data. When curve fitting from a transfer function, the modeling 
technique, no matter how precise, has to assume a kind of damp- 
ing behavior since modal damping is one of the parameters that has 
to be determined. This damping behavior is difficult to assess with 
precision and this is where we have an important source of error. 
The procedure developed here consists in taking a very high reso- 
lution drive point transfer function by zooming about each mode of 
interest, and then curve fitting by applying the nonlinear formula of 
the multimodal transfer function to a linear programming technique 
(Simplex) that minimizes the square of the residuals. 

6. The compatibility conditions used to connect structures normally es- 
tablish simple relations like x& = xb, but for some cases where we 
have a motor with lead screw, or just a motor M that is supported on 
one component A to drive another component B, we then use another 
relation which is x& + xm = ifl- This relation enforces the system 
kinematic relations and the force equilibrium as well. 

7. The use of the GPIB interface to establish communication between 
the HP analyzer and the PDP 11/44 computer really enhances the 
capabilities to do 3-dimensional analysis since it allows among other 
things use of an animated display of the computer predicted mode 
shapes. 

8. The eigenfunctions, determined experimentally for a number of dis- 
crete points along the ways on which the carriages slide , are subject 
to error and that adds bumps to the simulated ways. To smooth this 
we curve fit a beam-like eigenfunction to these data to get in addition 
to smoothness, a continuous function that permits a faster calctilation 
of the mode shape value at the point of connection (eliminating the 
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need for interpolation in real time). 

9. The use of the optical encoder motion that monitors the motor-rotor 
displacement during the robot operation to drive the equations of 
motion permits us to disregard the control system in the vibration 
analysis. After a computer-controlled experiment gives the encoder 
motion as well as the real vibration characteristics of the end effector, 
a computer simulation that utilizes the proposed mathematical model 
can be run using as input this encoder motion. 

10. A dynamic model for the system vibrations that normally contains as 
a forcing term the motor force can be manipulated to eliminate the 
unknown force in favor of the encoder motion (displacement, veloc- 
ity, and acceleration) and have this encoder motion drive the resulting 
system of coupled differential equations. This process is not restricted 
to time-invariant systems because it is done in the time domain, and 
it is simple to use for slowly-varying systems. The elimination of 
the unknown force which appears in every equation of the originally 
decoupled system is performed by combination of equations that in- 
volve multiplication by constants (from eigenvectors) and additions, 
plus the incorporation of an extra equation involving the encoder 
motion. 

11. The arrangement of computer programs to do the analysis is such that 
the simulation program receives the data in its simplest form, having 
already being processed by previous programs so that the real-time 
simulation can run efficiently. This fast program however is still 40 
times slower than the real-time process when running on a computer 
VAX/750 using 10 natural modes of vibration. 

12. The results of the mathematical model for Gxed configuration show 
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that the eigenvalues and eigenvectors correlate relatively well with 
the experimental values. The difference in natural frequencies for 
the range to 200 Hz is 6.89 % average, and this error is not so 
much affected by the convergence rate of the CMS which behaves 
like a regular Ritz analysis, but it is more affected by the errors in 
the modal parameter estimation, particularly in the scaling of the 
eigenvectors. 

13. The results of the mathematical model for the case of the moving 
robot show that the vibration follows very closely the encoder ac- 
celeration, and this in turn contains in its harmonics (among other 
harmonics introduced by the control system) those natural modes of 
the robot that are excited by the motion. A comparison of experi- 
mental and calculated vibrations show in general good correlation in 
magnitude and in frequency content. 

14. With regard to the preliminary work of the four-bar linkage, one anal- 
ysis that illustrates the nature of the position-dependent problems, or 
in other words, of changing boundary conditions is the rocker-beam 
analysis. Here we studied the behavior of the first natural frequency 
as the boundary conditions changed, having obtained the limiting 
cases of cantilever beam, pinned-pinned, clamped-pinned and free- 
pinned beams. 

15. Another result of the four-bar linkage problem, the linearization of 
the equations of motion by neglecting higher order terms leads to a 
parallel analysis in the moving 2-degree of freedom system on a beam 
which allowed extending the method of CMS to position-dependent 
problems. 

16. Finally, the dimensional analysis for the four-bar linkage illustrates 
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that dynamically similar systems behave in the same way, and there- 
fore, in the normalization of the mode shapes for the robot problem 
we have to distinguish the right member of the family which belongs 
to the data obtained in order to properly identify the system. 

5.2 Recommendations 

1. An improved technique for modal parameter identification should be 
pursued to more closely evaluate the modal parameters, particularly 
the orthonormal mode shape value since that determines most impor- 
tantly the precision of the final simulation results. Different damping 
patterns should be examined, and solution of the convergence prob- 
lems of multimodal curve fits should be tried. 

2. Another study should be. to incorporate the control system to the 
mathematical model in order to have a complete representation of the 
dynamic system. Again CMS may be used to connect the mechanical 
system to the control system through the motors. 

3. Another study which would make the model more precise would be to 
add corrections to first order to the differential equations by restoring 
some of the neglected terms. 

4. To make the program faster, one approach that could be tried is the 
asymptotic approximations to the solution of the differential equa- 
tions by curve-fitting first the appropriate data and then apply stan- 
dard techniques to approximately solve differential equations with 
slowly-varying coefficients. The solution would be expressed as a 
function of position and of the motion parameters. 

5. Another study should correct the mathematical model to account 
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Appendix A 



MODE SHAPE DATA 



A.l Determination of Modal Characteristics 

The vibration analysis to be performed on the cartesian robot requires the 
determination of the natural frequencies and mode shapes of a number 
of structures and components, some of them being rather complicated in 
shape. These parts were already described in Chapter 2 and include the 
robot structure, the carriages, and the lead screws. As mentioned before, 
the carriages were tested in a free condition (floating) by suspending them 
on soft foam-pads, and the impact test was applied to a number of test 
points. Also, the robot structure was tested while bolted to the floor by 
exciting it with an electrodynamic shaker which was suspended from the 
ceiling with a flexible string. Finally, the torsional modes of the lead screws 
were obtained analytically. These conditions are summarized in Table A.l. 
The random test using a shaker is more precise than the impact test 
when we need a calibrated measurement, and for the case of the robot 
structure attached to the floor, this was required in order to normalize the 
eigenvectors. For the Z-carriage we required more precision than for the 
other carriages since its modal characteristics tend to affect more the robot 
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vibrations. Also, its modes occur at higher frequencies when tested in free- 
free condition, and the impact test can not give good results above 800 Hz 
due to the roll-off of the impulse force spectrum. 

The test set-up for the impact procedure is shown in Figure A.l. The 
accelerometer is kept fixed at one convenient test point and the structure 
is marked at each of the test points. The instrumented hammer is used 
to impact on each of the test points while the analyzer collects a series of 
transfer functions. This process can be repeated for a different position of 
the accelerometer to get a new set of data which added to previous set could 
provide a more complete analysis; that is, the second set could show some 
of the modes more clearly than the first set. But strictly speaking, one set 
of tests is sufficient. The analyzer can then process all these data with an 
algorithm that determines the mode shape values at the test points. 

The test set-up for the random test is shown in Figure A. 2. Here the 
shaker is kept fixed to a convenient test point while the accelerometer has 
to be moved around for each new transfer function. 

The algorithm used by the analyzer (HP 5423A Structural Dynamics 
Analyzer) presents two options to the user, quadrature picking and rational 
fraction least squares, both for the single degree of freedom case. The 
quadrature picking or line-method is useful when damping is low and the 
modes are closely spaced since it only uses the transfer function value at the 
mode peak. The rational fraction least squares or band-method curve-fits 
a band of the transfer function that includes the mode being analyzed, and 
it is more precise when the modes are well separated. 

For the cases of very simple geometry, like the lead screws in torsion, 
we can analyze them by the theory of continuous-parameter models, par- 
ticularly when we have one or two uniform elements. Then the analysis is 
simple and accurate. 
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Figure A.l: Test set-up for impact procedure 
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Figure A. 2: Test set-up for random procedure 
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Figure A.3: Mode shapes no.l (top) and no. 2 (bottom) of robot structure 
as obtained by test 
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Figure A.4: Mode shapes no.3 (top) and no. 4 (bottom) of robot structure 
as obtained by test 
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Figure A. 5: Mode shapes no.5 (top) and no.6 (bottom) of robot structure 
as obtained by test 
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Figure A.6: Mode shapes no. 7 (top) and no. 8 (bottom) of robot structure 
as obtained by test 
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Figure A. 7: Mode shapes no. 9 (top) and no. 10 (bottom) of robot structure 
as obtained by test 
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Figure A. 8: Mode shapes no. 11 (top) and no. 12 (bottom) of robot structure 
as obtained by test 
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Table A.l: Conditions of experimental tests. 
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Figure A. 9: Mode shape no. 13 of robot structure as obtained by test 



APPENDIX A. MODE SHAPE DATA 



144 




50. 100. 

Frequency (Hs) 



150. 



200. 



Figure A. 10: Drive-point transfer function of robot structure 

A. 2 Mode Shape Data for Robot Structure 

The robot structure attached to the floor and without carriages gives the 
following 13 natural modes of vibration in the frequency range of to 200 
Hz (see Figures A.3- A.9). The drive-point transfer function which was 
used for normalization of the eigenvectors was obtained by zooming about 
each of the modes and then collecting all these transfer functions to make 
one single plot (see Figure A. 10). This figure shows a collection of nearly 
2000 points which give the transfer function absolute values with enough 
frequency resolution about each mode of interest. The curve fitting of these 
data was already discussed in section 2.3, and the curve fit is shown as a 
solid line in the same figure. From this curve fit we determine the modal 
parameters in a more precise way than in the standard modal analysis , 
and they are used then to normalize the eigenvectors. 
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A. 3 Mode Shape Data for Carriages 

All three carriages were tested separately in a free condition to get the 
elastic mode shapes shown in Figures A. 11- A. 20. The rigid body modes 
were calculated as explained in Chapter 2 . 

The normalization was carried out in two different ways according to 
the modal test used. For the Z-carriage we used the random test and from 
it we got the drive-point transfer function shown in Figure A. 21 which is a 
collection of zoom transfer functions about each mode. We then curve-fit 
the modal parameters as explained in section 2.3. 

The X-carriage and Y-carriage were tested by the impact procedure (un- 
calibrated) and, then to normalize we used the mass distribution according 
to a formula of Chapter 2. 

A. 4 Mode Shapes for Lead Screws 

These mode shapes in torsion were obtained analytically by modeling each 
lead screw as a uniform shaft with an inertia (motor rotor) attached to one 
end, as shown in Figure A.22. The natural frequencies for this system are 
given by the transcendental equation 

Joj [G 



GI P 



y|=-tan(«Z^) (AA7) 

and the eigenfunctions are given by 

a(x) = Cicos(ojxJ— ) + C2sin(ujxJ— ) (A. 48) 

To normalize we use orthogonality of the normal modes to get 

f L P I t> al{x)dx + Ja 2 n {L) = 1 (A49) 
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Figure A. 11: Mode shapes no.l (top) and no.2 (bottom) of X-carriage as 
obtained by test 
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Figure A. 12: Mode shapes no.3 (top) and no. 4 (bottom) of X-carriage as 
obtained by test 
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Figure A. 13: Mode shapes no.l (top) and no. 2 (bottom) of Y-carriage as 
obtained by test 
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Figure A. 14: Mode shapes no.3 (top) and no. 4 (bottom) of Y-carriage as 
obtained by test 
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Figure A. 15: Mode shapes no. 5 (top) and no. 6 (bottom) of Y-carriage as 
obtained by test 
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Figure A. 16: Mode shapes no. 7 (top) and no. 8 (bottom) of Y-carriage as 
obtained by test 
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Figure A. 17: Mode shapes no. 9 (top) and no. 10 (bottom) of Y- carriage as 
obtained by test 
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Figure A. 18: Mode shapes no.l (top) and no. 2 (bottom) of Z-carriage as 
obtained by test 
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Figure A. 19: Mode shapes no. 3 (top) and no. 4 (bottom) of Z-carriage as 
obtained by test 
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Figure A. 20: Mode shape no. 5 of Z-carriage as obtained by test 
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Figure A. 21: Drive-point transfer function of Z-carriage 
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I p Polar moment of area of shaft 
J Moment of inertia of rotor 

Figure A. 22: Model of motor-rotor and lead-screw 

The results for both lead screws were given in Table 2.3. The Z-cairiage 
motor MZ with its short shaft and pinion is essentially rigid, so we consid- 
ered only its rigid body mode in torsion. 



A. 5 Modal Matrices for Data Input to Math- 
ematical Model 

The modal matrices to be used in the mathematical model are obtained 
from the mode shapes just shown. To get these data files into the computer, 
we transfer data from the analyzer through a GPIB interface by means of 
a computer program shown in Table A. 2. 

The modal matrices are formed from these data files by properly scal- 
ing them in order to normalize the eigenvectors. If the correct drive-point 
mode-shape value has been obtained by curve-fitting the drive-point trans- 
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Table A. 2: Computer programs to transmit data, back and forth, between 
a pdpll/44 computer and an HP analyzer. 

C TRAMOD.FTN 

(i************************* 

C 

PROGRAM TO TRANSFER HP5423A MODAL DATA TO THE 

C PDP-11/44 VIA THE NATIONAL INSTRUMENTS GPID11-1 INTERFACE. 

C 

C DATA IS TRANSFERRED IN ASCII DY A 501,4.2 SAVE COMMAND 
C 

C THE OUTPUT FILE IS FORMATTED WITH 1 VALUE PER RECORD (E15.7) 

C NOTE: THE IB: HANDLER USES RSX CHANNEL 1 
C 

C SEE CHAPTER 3 OF VOL. 3 "USING HP-IB" (5423A) 

£**»******«»**♦**♦******»«*»*» 

C 

BYTE BDATA(630), BSAVE(IO), FILNAM(20),CSAVE(96) 

REAL MDATA(999) 

INTEGER*4 ASAVE(3) 

C 

C 501,4,2SA IS THE SAVE TO CONTROLLER COMMAND 

C 10 IS THE DECIMAL VALUE OF THE LF TERMINATOR. 

C 

DATA BSAVE / '570717,747,72','S','AM0/ 

DATAIWR,IRE,ICL.ITRJRM,n,0,ffO,ICO,B?C i rDE,IFI 

1 /0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10/ 

DATA fflP, ITA, ILI /"4, "104, "44/ 

C 

C TELL THE OPERATOR TO PREPARE THE INTERFACE 
C 

TYPE *,' TRAMOD- HP5423A to PDP-11/44 data transfer program.' 

TYPE *,' ' 

TYPE VI. Connect the HP-IB cable from the PDP-lltotheHP.' 

TYPE *,'2. Sot the HP back panel switch to "ADDRESSABLEONLYV 

5 TYPE V3. Get measurement data on the active trace of the HP.' 

TYPE *,'4. Enter to stop or 1 to continue :' ■ 

ACCEPT *, J 

IF(J .EQ. 0) GO TO 1000 

TYPE.*,' 5. ENTER MODE, 1ST DOF, AND LAST DOF: ONE AT A TIME' 

ACCEPT *, ASAVE(l) 

ACCEPT *, ASAVE(2) 

ACCEPT *, ASAVE(3) 

NWRDS=ASAVE(3)-ASAVE(2)+1 
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Table A. 2: continued 

C 

C FIRST CLEAR THE BUS 

C 

10 TYPE *,' Initialising ...' 

J = IDUP(ICL,-1) 

IF( J .NE. 1)TYPE Y BUS CLEAR FAILED, ERROR=\J 

C 

C MAKE SURE THAT ANALYZER IS DEFINED AND ON THE BUS. 

C 

J = IBUP( IDE, IHP, ITA, ILI, 0, 0, 0, 0) 

IF( J .NE. 1)TYPE *,' HP DEFINE FAILED, ERROR=',J 

J = IBUP( ICL, IHP) 

IF( J .EQ. l)GO TO 20 

C 

15 TYPE V TRY AGAIN ? 0=NO, STOP; 1=YES :' 

ACCEPT *,J 

IF(J .EQ. 0) STOP 

GO TO 10 

C 

C ALL SET, SEND 501,4,2 SAVE COMMAND 

C 

20 ILEN = 10 - 

J = IBUP( IWR, IHP, BSAVE(l), ILEN) 

IF(J .EQ. ILEN)GO TO 30 

C 

C ERROR 

C 

25 TYPE *,' SAVE COMMAND FAILED, REC=',J,' SENT=',ILEN 

C 

C SERIAL POLL TO GET SDA SRQ 

C 

J = IBUP( IPO, IHP) 

TYPE *,' ANALYZER SRQ = ',J 

GOTO 15 

C 

C START READING HEADER (3 VARIABLES) 

C 

C THIS POLL IS NEEDED WHEN CONTROLLER SENDS SAVE COMMAND 

C THE PROGRAM WAITS UNTtt THE ANALYZER IS READY (SRQ=96) 

C 

30 ISAV = 91 

DO 35 I = 1, 2000 

J = IBUP(H > 0, IHP) 

B?(J .EQ. ISAV)GO TO 38 
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Table A. 2: continued 



35 CONTINUE 

TypE * , TIME0 UT ON SAVSRQ, HP SRQ =\J 

GO TO 25 

C 

38 ILEN = 32 

IOF = 1 

ILAN=96 

ENCODE(ILAN,891,CSAVE(l))ASAVE(l),ASAVE(2),ASAVE(3) 

891 FORMAT(3I5) 

J=IBUP(IWR,IHP,CSAVE(1),ILAN) 

C 

C NOW READ IN NWRDS OF (REAL) DATA TO MDATA 

C USE BDATA( 17-33) AS A BUFFER 

C 

999 FORMAT(E14.4) 

ILEN = 32 

DO 50 I = 1, NWRDS 

J = IBUP( IRE. IHP, BDATA(17), ILEN) 

DECODE( ILEN-2, 999, BDATA(17) )MDATA(I) 

TYPE *,MDATA(I) 

50 CONTINUE 

C 

C STORE ON DISK FILE IN ASCII (2 VALUES / RECORD) 

C 

60 TYPE 996 

996 FORMATC Enter filename for storage :',$) 
ACCEPT 995,FILNAM 

995 FORMAT(20A1) 

FDwNAM(20) = 

C 

OPEN( UNIT=2, NAME=FILNAM, TYPE='NEW', ERR=60) 

C 

WRITE(2,994) (MDATA(I),I=1,NWRDS) 

994 FORMAT(E15.7) 

CLOSE(UNIT=2) 

TYPE 997,FttNAM 

997 F0RMAT(1X/Ascii data and header stored in "\20A1,"") 
GO TO 5 

1000 STOP 
END 
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Table A. 2: continued 

rnt**^**t******* ************************* *********************** 

C TRARES2.FTN 

C PROGRAM TO TRANSFER HP5423A RESIDUE DATA FROM THE 

C PDP-11/44 VIA THE NATIONAL INSTRUMENTS GPID11-1 INTERFACE. 

C 

C DATA IS TRANSFERRED IN ASCII BY A 504,1 RECALL COMMAND 

C 

C THE INPUT FILE SHOULD BE FORMATTED WITH 1 VALUE PER RECORD (E15.7) 

C 

C THE PROGRAM TAKES THE 40 VARIABLE HEADER FROM THE END OF THE FILE 

C AND SENDS IT FIRST. 

C 

C NOTE: THE IB: HANDLER USES RSX LOGICAL UNIT 1 

C SEE THE NATIONAL INSTRUMENTS "GPIB11 SOFTWARE REFERENCE MANUAL" 

C FOR INFORMATION ON USING THE BUS FROM FORTRAN OR ASSEMBLER 

C 

C SEE CHAPTER 3 OF VOL. 3 "USING HP-IB" (5423A) 

C FOR THE SAVE/RECALL DATA AND HEADER FORMATS 

C 

p*» ************************************************* ***************** 

C 

BYTE BDATA(630), BRCL(IO), FttNAM(20),CSAVE(96) 

REAL MDATA(700) 

INTEGERS ASAVE(3) 

C 

C 504,1RA IS THE RECALL FROM CONTROLLER COMMAND 

C 10 IS THE DECBvIAL VALUE OF THE LF TERMINATOR. 

C 

DATA BRCL / '5\'0\T,7, , 47,71\ , R7A',10/ 

C 

C mSTRUMENT BUS UTDLITY PROGRAM (D3UP) FUNCTION CODES 

C 

DATArWR.IRE,ICL,ITR,IRM,n.O,n , 0,ICOJPC4DE,IFI 

1 /0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10/ 

C 

DATA B3P, ITA, HJ /"4, "104, "44/ 

C 

C TELL THE OPERATOR TO PREPARE THE INTERFACE 

C 

TYPE *,' TRARES2- HP5423A from PDP-11/44 data transferprogram.' 

TYPE V ' 

TYPE *,'l. Connect the HP-B3 cable from the PDP-11 to theHP.' 
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Table A. 2: continued 

TYPE \'2. Sot the HP back panel switch to "ADDRESSABLEONLYV 

10 TYPE *,'3. Enter to stop or 1 to continue :' 

ACCEPT *,J 

IF(J .EQ. 0) GO TO 1000 

TYPE *, ' 4. ENTER MODE, 1ST MEAS AND LAST MEAS [MAX MEAS+4IF' 

TYPE VMODE.NE.0], [(MAX MEAS+1)*4 IF MODE.EQ.0], ONE BY ONE' 

ACCEPT *, ASAVE(l) 

ACCEPT *, ASAVE(2) 

ACCEPT *, ASAVE(3) 

NWRDS=ASAVE(3)-ASAVE(2)+1 

C 

C FIRST CLEAR THE BUS 

C 

20 TYPE *,' Initialising ...' 

J = IBUP(ICL,-1) 

IF( J .NE. 1)TYPE *,' BUS CLEAR FAILED, ERROR=\J 

C 

C MAKE SURE THAT ANALYZER IS DEFINED AND ON THE BUS. 

C 

J = IBUP( DDE, IHP, ITA, ILL. 0, 0, 0, 0) 

D7( J .NE. 1)TYPE *,' HP DEFINE FAILED, ERROR=\J 

J = IBUP( ICL, D3P) 

W{ 3 J2Q. l)GO TO 40 

C 

30 TYPE V TRY AGAIN ? 0=NO, STOP; 1=YES :' 

ACCEPT *,J 

D?(J .EQ. 0)GO TO 1000 

GO TO 20 

C 

C GET NAME OF ASCH FILE (1 VALUE / RECORD) 

C 

40 TYPE 999 

999 FORMATp Enter filename for transfer (0=quit) :',$) 

ACCEPT 998,Fn.NAM 

998 FORMAT(20A1) 

IT = FH.NAM(1) 

IF(IT .EQ. 48) GO TO 1000 

FILNAM(20) = 

C 

C GET NUMBER OF DATA VALUES TO SEND 

C 

50 TYPE 997 
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Table A. 2: continued 



997 FOB.MATC Enter number of data values to be sent -',/, 

1' 512 for trans or time, 250 for coh :',$) 

ACCEPT *, NVAL 

C 

C READ DATA AND HEADER FROM FILE 

C 

D DO 45 1=1,630 

D 45 BDATA(I) = 

OPEN( UNIT=2, NAME=FILNAM, TYPE='OLD', ERR=40) 

C 

READ(2,996) (MDATA(I),I=1,NVAL) 

996 FORMAT(E15.7) 

CLOSE(UNIT=2) 

C 

C CHECK THAT NWRDS IN HEADER (VARIABLE 3) = 2 * NVAL 

C THE ENCODE ROUTINE TRANSLATES FROM DEC INTERNAL FORM TO ASCH 

C THE DECODE ROUTINE TRANSLATES FROM ASCII TO DEC INTERNAL FORM 

C THE LAST TWO CHARACTERS ARE CR,LF WHICH ARE NOT DECODED. 

C 

C ILEN = 16 

C IOF = ILEN * 2 + 1 

C DECODE(ILEN-2, 995, BDATA(IOF) ) BUF 

C995 FORMAT(E14.4) 

C NWRDS = BUF 

C 1F(NWRDS .EQ. 2*NVAL) GO TO 60 

C TYPE *,' ** ERROR, NVAL DISAGREES WITH HEADER' 

C GO TO 60 

C 

C ALL SET, SEND 501,4,1 RECALL COMMAND 

C 

60 HEN = 10 

TYPE *,' Sending 501,1 recall command...' 

J = B3UP( 1WR, HIP, BRCL(l), ILEN) 

IF(J .EQ. E-EN)GO TO 80 

C 

C ERROR 

C 

70 TYPE *,' RECALL COMMAND FAILED, REC=\J,' SENT=',H J EN 

C 

C SERIAL POLL TO GET SDA SRQ 

C 

j = roup( ipo, mp) 
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TYPE *,' ANALYZER SRQ = \J 

GO TO 30 

C 

C THIS POLL IS NEEDED WHEN CONTROLLER SENDS RECALL COMMAND 

C THE PROGRAM WAITS UNTIL THE ANALYZER IS READY (SRQ=121) 

C 

80 IRCL = 121 

DO 90 I = 1, 2000 

j = mur(iPO, rap) 

IF(J .EQ. IRCL)GO TO 100 

90 CONTINUE 

TYPE *,' TIMEOUT ON RECALL SRQ, HP SRQ =',J 

GO TO 30 

C 

C START SENDING HEADER (3 VARIABLES) 

C 

100 CONTINUE 

TYPE *,' Sending measurement header...' 

ILEN=32 

IOF=l 

ILAN=96 

ENCODE(ILAN,891,CSAVE(l))ASAVE(l),ASAVE(2),ASAVE(3) 

891 FORMAT(3I5) 

J=IBUP(IWRJHP,CSAVE(1),ILAN) 

C ILEN = 16 

C IOF = 1 

C J = IBUP( rWR, fflP, BDATA(IOF), ILEN) 

C IF(J .NE. ILAN)TYPE *,' FIRST VARIABLE TRANSMIT FAILED, J=',J 

C 

C SEND THE REST OF THE HEADER (VARS 2-37) 

C 

C DO 110 I = 2, 37 

C IOF = IOF + ILEN 

C TYPE 333,(BDATA(I0F+J-1),J=1JLEN) 

333 F0RMAT(22A1) 

C J = B3UP( IWR, IHP, BDATA(IOF), ILEN) 

110 CONTINUE 

C 

C LAST THREE VARS ARE ASCH 

C 

C IOF = IOF + ILEN 

C ILEN = 8 
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C TYPE 333,(BDATA(IOF+J-l),J=l,ILEN) 

C J = IBUP( IWR, IHP, BDATA(IOF), ILEN) IINPUT TRANSDUCER MODEL 

C IOF = IOF + ILEN 

C TYPE 333,(BDATA(IOF+J-l),J=l,ILEN) 

C J = IBUP( IWR, IHP, BDATA(IOF), ILEN) !OUTPUT TRANSD. 

C IOF = IOF + ILEN 

C ILEti = 22 

C TYPE 333,(BDATA(IOF+J-l),J=l,HJ3N) 

C J = IBUP( IWR, IHP, BDATA(IOF), DLEN) !MEAS. TITLE 20 CHARS. 

C TYPE 994,(BDATA(IOF+I-l)a=1.20),NVAL 

994 F0RMAT(1X,' Title= "\20A1,"",/,' Sending',14,' datavaliMs.') 

C 

C NOW SEND OUT NVAL OF (REAL) DATA TO THE ANALYZER 

C USE BDATA(17-32) AS A BUFFER 

C FINISH EACH VARIABLE WITH CR,LF FOR TOTAL OF 16 BYTES / VARIABLE 

C 

BDATA(31) = 13 

BDATA(32) = 10 

ILEN = 16 

993 FORMAT(E14.7) 

DO 120 I = 1, NVAL 

ENCODE( n.EN-2, 993, BDATA(17) )MDATA(I) 

TYPE *,MDATA(I) 

J = D3UP( IWR, fflP, BDATA(17), HJEN) 

120 CONTBVUE 

TYPE *,' Operation finished.' 

GO TO 10 

1000 CONTBWE 

STOP 

END 
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fer function, then we scale the eigenvector by pivoting on the current drive- 
point mode-shape value. If the normalization procedure is based on the 
mass distribution of the test points of the component, then we determine 
C, the scaling factor, directly. 

Each robot component was tested on a number of points which include 
the points of connection to other components. These test points are marked 
all over the contour of each component in order to see how the mode shapes 
are. For data input to the mathematical model, we only need the points of 
connection. 

A. 6 Mode Shape Data for Complete Robot, 
Experimental and Predicted Results 

The robot was tested for the configuration shown in Figure 4.1, that ia, 
with the X-carriage and Y-carriage in the middle and the Z-carriage down. 
The mode shapes in the frequency bandwidth of to 200 Hz, as obtained by 
experimental test using a shaker, are shown in the top part of Figures A.23- 
A.39. 

For this same configuration, we have run the mathematical model to 
predict the mode shapes of the robot; the synthesized mode shapes are 
shown in the bottom part of Figures A. 23- A. 39. A comparison of natural 
frequencies as well as modal damping factors can be seen in Table A.3. 
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Figure A. 23: mode shape no.l of complete robot as obtained by test(top) 
and by synthesis(bottom) 
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Figure A. 24: mode shape no. 2 of complete robot as obtained by test(top) 
and by synthesis(bottom) 
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Figure A. 25: mode shape no. 3 of complete robot as obtained by test(top) 
and by synthesis (bottom) 
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Figure A.26: mode shape no.4 of complete robot as obtained by test(top) 
and by synthesis(bottom) 
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Figure A. 27: mode shape no. 5 of complete robot as obtained by test(top) 
and by synthesis(bottoin) 
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Figure A. 28: mode shape no.6 of complete robot as obtained by test(top) 
and by synthesis(bottom) 
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Figure A. 29: mode shape no. 7 of complete robot as obtained by test(top) 
and by synthesis(bottom) 
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Figure A. 30: mode shape no. 8 of complete robot as obtained by test (top) 
and by synthesis(bottom) 
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Figure A. 31: mode shape no.9 of complete robot as obtained by test(top) 
and by synthesis(bottom) 
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Figure A.32: mode shape no. 10 of complete robot as obtained by test(top) 
and by synthesis(bottom) 
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Figure A.33: mode shape no. 11 of complete robot as obtained by test(top) 
and by synthesis (bottom) 
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Figure A. 34: mode shape no. 12 of complete robot as obtained by test(top) 
and by synthesis(bottom) 
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Figure A.35: mode shape no.13 of complete robot as obtained by tcst(top) 
and by synthesis(bottom) 
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Figure A.36: mode shape no.14 of complete robot as obtained by test(top) 
and by synthesis (bottom) 
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Figure A. 37: mode shape no.15 of complete robot as obtained by test(top) 
and by synthesis(bottom) 
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Figure A. 38: mode shape no. 16 of complete robot as obtained by test(top) 
and by synthesis (bottom) 
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Figure A. 39: mode shape no. 17 of complete robot as obtained by test(top) 
and by synthesis(bottom) 



APPENDIX A. MODE SHAPE DATA 



183 



Mode 
No. 


Meas 
/» (Hz) 


ured 
?« (%) 


Predicted 
/„ (Hz) [ fc (%) 


1 


12.36 


1.98 


11.93 


1.45 


2 


18.85 


1.65 


19.17 


1.45 


3 


36.36 


1.42 


36.87 


0.78 


4 


42.83 


9.52 


46.34 


3.30 


5 


47.83 


4.90 


49.53 


2.54 


6 


60.51 


3.24 


61.38 


1.45 


7 


70.72 


0.89 


75.06 


1.82 


8 


81.47 


1.95 


96.47 


0.63 


9 


86.67 


1.74 


101.75 


1.41 


10 


97.70 


1.17 


123.39 


1.30 


11 


112.52 


7.39 


126.84 


1.06 


12 


126.27 


0.10 


128.72 


0.75 


13 


139.43 


2.62 


141.05 


1.16 


14 


149.85 


1.00 


145.68 


2.86 


15 


152.34 


0.66 


148.23 


2.76 


16 


178.12 


1.79 


169.39 


1.04 


17 


184.37 


0.01 


179.06 


1.11 



Table A. 3: Comparison of natural frequencies and modal damping factors 



Appendix B 

VIBRATION OF MOVING 
ROBOT 



This appendix describes the method of analysis employed in the study of 
robot vibrations during motion,. method which has been used in Chapter 3 
to do the simulations. 

B.l Methods of Analysis 

After briefly mentioning the case of prescribed rigid body motions which 
provides a simple approach to the problem, we then discuss in some detail 
the procedures that can be used for the case when the motor-rotor vibration 
is specified. 

B.l.l Prescribed Rigid Body Motion 

As was discussed in Chapter 3, one of the simplest approaches to the prob- 
lem of predicting the robot vibrations caused by its motion is to prescribe 
how a particular substructure is moving (as a rigid body) regardless of the 
input required to cause that motion, thus leaviug out the effects of the robot 

184 
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control system. In order for this procedure to be simple we can assume for 
a moment that the substructures are rigid and calculate the forces exerted 
at the connection points; forces which would produce that specified motion 
of the center of gravity of the stibstructure under the presence of friction 
and gravity. This procedure, although a rough approximation leads to a 
relatively fast computer model. 

B.1.2 Prescribed Motor- Rotor Vibration 

The motor-rotor motion can be easily determined experimentally from en- 
coder data, and this motion, including vibration, may be used as input to 
the mathematical model. We drive the mathematical model using this pre- 
scribed motor-rotor vibration while ignoring the forcing function exerted 
by the amplifier- mo tor combination. The procedure developed for this pur- 
pose had to be restricted to the available data, that is, to the use of modal 
parameters only. Two suitable techniques for this analysis are described in 
what follows. 

Frequency-Domain Technique 

There is a very powerful technique for this application which works in the 
frequency domain through the use of Fourier transforms, unfortunately it 
is only valid for time-invariant systems. This can be briefly described as 
follows: Consider a flexible system as shown in Figure B.l where a force 
f p (t) is applied at the point P producing a vibration x p (t). Suppose we 
know x p (t) and the transfer function 



*-<"> = tw (BM) 
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Flexible system 



Figure B.l: Flexible system 
which may be synthesized from the modal parameters as follows: 

Then, by taking the Fourier transform of x p (t) to get X p (uj) we can calculate 
F p (<jj), and after taking the inverse Fourier transform we can get f p (t). With 
f p (t) known, we can then integrate the equations of motion 

{p} + [2f nWn ]{p} + [ w »]{p} = [*]«{/> (B.52) 

to get the vibration at any point of the structure by means of {xq} = 

Another procedure which may be faster to apply consists of the follow- 
ing. Assume that we also know Hpq(uj) as indicated below, and we wish 
to calculate ig(t). 

HpqH = £ ,/Vr V- (B.53) 
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Then 

t v\ u} ) 
Solving for F p (u>) in equation B.50 and B.54 we can then equate both 

expressions to get: 

Then we calculate x<y(£) from equation B.55 by inverse transforming Xq(u>) 
, all without explicitly calculating the force fp(t). 

One source of error in this process may come from the use of the FFT 
which assumes that the function xp{t) is periodic, but in this case of robot 
motions there is a start point. This start condition is normally taken into 
account by extending the time history back from time zero, that is: 



x P (t) = 



for t < 

(B.56) 
x P {t) for t>0 



and giving enough span to get good approximation. 

Even though this technique was not applied to the study of the robot 
vibrations, it could be used to get approximate results if we considered 
average values for the transfer functions involved. 

Time-Domain technique 

A time-domain technique which can be of more general application is the 
following: 

Consider the same system of Figure B.l where we ultimately want to 
determine the vibration x(t) at any point by knowing xp(t) but not fp(t). 
The system can be described by the truncated system of uncoupled equa- 
tions B.52. Since there is only one external force fp(t) acting on the system, 
we can express the modal force as 

mn = w'{/p> (3.57) 
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where {a} is formed by the modal matrix evaluated at point P in the 
direction of the force. 

In addition to the n equations of B.52, we have another condition to be 
satisfied: 

M = \M{P) = {«}'{?} (5-58) 

So, we have n + 1 equations and n + 1 unknowns, those are {p} and fp. 

Assuming that a, ^ for all i we can eliminate from equation B.52 fp 
by combining the n equations to get n — 1 equations as follows: 

• multiply 1st equation by £ : ±pi + ^-fiWxPi + ^-w?Pi = fp- 

• multiply i + 1th equation by -~ : 

Subtract to get the ith new equation without the term in fp. Then we 
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get n — 1 equations for the n variables p, 
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Now we can apply equation B.58 to add another equation (or row) in 
such a way that the mass matrix remains non-singular, that is, by differ- 
entiating equation B.58 twice with respect to time: 



{x~p} = [4>p]{p} = { a y{p} 



(5.60) 



so, the nth equation of B.59 after the addition would contain {a} in the mass 
matrix and the corresponding forcing function would be {xjr>}« Therefore, 
starting with a coupled system of equations given by equation B.52 and 
equation B.58, we end up with a coupled system of equations which is 
driven by the prescribed motion xp(t), and which can be integrated directly 
to determine {p} and ultimately {x}. 
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Substructure A Substructure B 

Figure B.2: Two substructures A,B 

One important observation regarding the application of this method is 
that the differentiation of xp with respect to time assumes that fa = 0. 
This is strictly correct for a time-invariant system, but only approximately 
valid for a slowly-varying system. For the general case we would have to 
include terms in fa and fa in the equations of motion. 

APPLICATION TO A SIMPLE CASE 

In order to evaluate these methods, we consider a simple case which can 
illustrate with detail what we have applied to the study of the cartesian 
robot. To make it more similar to the robot case, we include a connection 
of two substructures. 

The problem is defined as follows. Consider the two substructures, A 
and B, as shown in Figure B.2. Given xi(t) of the assembled system, find 
x 2 (i) of the same assembled system for two different cases: 



APPENDIX B. VIBRATION OF MOVING ROBOT 



191 



• Impulse response. 



Step response. 



Synthesis of Modal Data for Assembled System 
Substructure A has the following modal parameters: 



\Ua\ = 
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Wl 



1 = ^,[^b] = 4=,{/b} = o 



m 



■v/m' 



The uncoupled equations before assembly are given by 



{Pa} + [* A }{Pa} = [+A) T {fA} = { 
{Pfl} + [«Sl{pfl} = [<^]{/b} = 



\/2m 
. y/2m 



{0} 



(5.61) 



(5.62) 



(5.63) 



The compatibility condition requires that 



{x 2 } = {a*}, or [4> Ai ]{pa} = l^HPfl} (5.64) 



Choosing p# as the generalized coordinate to be eliminated, we get 



{PB} = [<t>B}- 1 l4>A 1 ]{pA} = [± ~^}{Pa} 



(5.65) 
Then, the transformation matrix T as defined by formula 2.7 becomes 



[T} = 



1 







' 





1 


, for the case < 


9i = PAi 


1 

v'2 


1 

V2 . 




92 = PAj 



(5.66) 
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The assembled system is obtained by 

[m] = [t\ t {t\, [k] = [r] vim in = mwm 

Substituting the previously obtained values, we get: 



_i 
2 

1 3 

2 2 



9i 
92 



> + 



2 2 



> = < 



\Z2fi 

_£_ 

1 v/2m 



. This assembled system has the eigenvalues 



and the eigenvectors 



[0] = 



011 012 
021 022 



(B.67) 



(B.68) 



'4 m _ V jt» 2 



(B.70) 



This can be verified to be exactly correct for the resulting system shown 
in Figure B.3. 



Analysis by Frequency-Domain Technique 

In order to get the response at xi given the response at ii, we can apply 
the frequency-domain technique described above. First we decouple the 
system of equations given by B.68 to get 



1 


- 


' 


u 


■«? 





< 




. = . 


1 




r 2 


1 





«ij 




r-i 





(0n+02i);fc 

. (012 + 022) -jt . 



. The transfer function at point x\ is 

0?i 



Hu(u) = 



u\ — u> 2 



+ 



2 



12 



w|-w 2 



(B.71) 



(B.72) 



So, by applying equation B.72 we can obtain Fi(u) and from there fi(t). 
Integrating equations B.68 we can determine Xi(t) as shown below, or we 
can use Hu in combination with Hu in formula B.55 to get ^(w) directly, 
and from it x 2 (t). 
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Figure B.3: Assembled A-B 
Analysis by Time-Domain Technique 
The prescribed vibration Xi(t) is related to the variables q by 

9i 



'i-WM-ijb ^} 



92 



(B.73) 



. First we eliminate j\ from equation B.68 by subtracting the first row 
from the 2nd to get 



-2 2 



4i 



92 



> + 



-k B 2k A + k B 



9i 

92 



= 



(fl.74) 



. Differentiating equation B.73 twice with respect to time, and adding that 
equation to B.74 we get 



-2 2 

l i 

v2m v'2m 



( \ 



9i 

92 



> + 



-k B 2k A + k D 


< 


Ul 1 

r = i 







92 J 



(B.75) 



Xl 
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These equations B.75 can be integrated directly to get qi(t) and (72 (0 
to then get 

i 2 (0 = [<fo]{q} = 4>2i1i + 4>2*q2 (B.76) 

An alternative way to arrive at the differential equation driven by the 
prescribed motion xi consists of eliminating one of the qi , say q 2 from 
equation B.73, and substitute in equation B.75 to get a single equation 

91 + -^r qi = V~ Xl + 4^r {2kA + ko)Xx {BJ7) 

Results 

Two kinds of inputs or prescribed motions at xi can be employed now. The 
impulse response is given by 

A2 ±2 

x 2 (t) = 22! sin (u 2 t) + ^i sin (u^i) (5.78) 

U) 2 Wj 

and the step response by 

±2 i2 

x 2 {t) = 1^.(1 - cos (ujit)) + r» (l - cos (w 2 0) (£.79) 



The corresponding responses at x 2 as given by the first method are 
shown in Figure B.4 and Figure B.5, where the calculated response is the 
top plot and the expected response 22(0 is shown for comparison at the 
bottom. Both results, the calculated and the expected, are very similar. 
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Figure B.4: Response at x% to impulse- response input at x\ t calculated 
(top) and expected (bottom) 
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Figure B.5: Response at i 2 to step-response input at ii, calculated (top) 
and expected (bottom) 



Appendix C 

PRELIMINARY STUDY OF 

VIBRATIONS OF 

A GEOMETRICALLY 

NONLINEAR SYSTEM 



C.l Preliminary Vibration Analysis 

Before we started doing the analysis of the cartesian robot, we were involved 
in the study of the vibration characteristics of a four-bar linkage with elastic 
input and output shafts which presented some common characteristics with 
the cartesian robot, and from whose analysis we could get insights into the 
problem of position-dependent systems. The experience obtained in the 
preliminary study facilitated the development of a procedure to analyze 
the robot vibrations, and that is why we dedicate this Appendix to its 
presentation. 

The system consists of a crank- rocker mechanism with flexible input 
and output shafts (sec Figure C.l). The input shaft transmits power from 



197 



AI'I'VNDIX C. PRELIMINARY STUDY OF VIllllATIONS 

1TI2 , l_2 ,Kg 



198 



m 3' L 3 
E 3» r 33. A 3 







Figure C.l: Continuous parameter model of the system 

an ideal motor, represented by its rotor inertia [ , to the crank. The input 
shaft stiffness is normally high as is the first natural torsional frequency 
of the shaft considered in isolation. Also the shaft inertia is usually quite 
small as compared to that of the rotor. Therefore, the shaft may be modeled 
as a linear torsional spring and viscous damper with spring and damping 
coefficients K\ and C\. 

For this study the crank and connecting rod are assumed to be short 
and their natural frequencies in bending, considered as isolated beams are 
assumed high when compared with the length and bending natural fre- 
quency of the rocker. Therefore the input link is considered to be a rigid 
link of mass mi and inertia Ii and the coupler link is modeled as a link of 
mass m 2 inertia I 2ct and axial stiffness K 2 . The rocker is a long uniform 
member; it has mass m 3 , and elasticity in bending only. The output shaft 
is rigidly attached to the rocker, and is perpendicular to the rocker's plane 
of motion. It is a uniform circular shaft of length L M , diameter d„, polar 
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moment of inertia I pa , density p„ ,and shear modulus G a . At the end of the 
shaft there is a system load represented as an inertia 74. 

Before proceeding to derive the equations, we should point out the sim- 
ilarities and differences with the robot problem: 

• For a given configuration of these systems, the modal characteristics 
can be calculated by one of several methods, and these modal param- 
eters may be different from those corresponding to another configura- 
tion, that is, by changing relative positions the boundary conditions 
of some elements may change, or the relative inertias and stiffnesses 
may vary. This is true for both systems, and as we will see later, the 
method to obtain the modal parameters has to be chosen according 
to the system complexity, and the more complex the system, the sim- 
pler the method ought to be in order to be able to solve the resulting 
differential equations. 

• The four-bar linkage problem consists of finding steady-state solutions 
and the cyclic variation of the system coefficients is fast as compared 
to the vibration which takes place superposed on the rigid body mo- 
tion. On the other hand, the cartesian robot problem concerns finding 
transient vibrations which occur during motion of the carriages; tiffs 
system can be modeled as a slowly varying one, that is, the system 
coefficients vary slowly as compared to the vibrations. In both cases, 
the solution involves integration of a system of differential equations 
with variable coefficients. 

C.2 Introduction 

Many authors have studied the vibration of four-bar linkages with one or 
more clastic links. In 1972, Erdman, Sandor and Oakberg [21] reported 
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using Finite Element methods to study an clastic linkage. In their analysis 
they considered effects of input speed fluctuations. In 1973, Imam, Sandor 
and Kramer [42] applied quasi-static structural techniques to this prob- 
lem; they included in their analysis the rate of change of eigenvalues and 
eigenvectors to reduce the computer time. In 1973, Sadler and Sandor [68] 
analyzed a linkage with rigid crank and elastic coupler links using Euler- 
Bernoulli theory for beams; to solve the equations they used the difference 
approximation and Taylor expansions. The same authors in 1974 [69] ana- 
lyzed a crank- rocker mechanism with rotational inertia in the output. They 
modeled the crank link as a cantilever beam and the coupler and rocker 
links as simply supported beams. Sutherland [74] , in 1976 assumed sinu- 
soidal mode shapes for the elastic links in order to model their behavior; 
he applied both analytical and numerical methods to solve the equations. 
For all of these analyses, the input link is assumed to be driven directly by 
the motor; for many cases, input link angular velocity is considered con- 
stant and there are no loads assumed to be acting on the system except 
the inertia loads of the links. In fact however, most linkages in use are 
driven through input shafts which are not rigid and motions of interest 
are often those of inertias being driven by shafts connected to the output 
or rocker link. Examples are found in industrial sewing machines, paper 
making machines, packaging machines, and others. Often the frequencies 
of the first torsional resonances of the input and output shafts are near 
the bending frequencies of the links. As a result, the resonances of the 
system are coupled and cannot be evaluated correctly unless the character- 
istics of the input and output shafts are taken into consideration. In 1978, 
Sanders and Tesar [70] performed an analytical and experimental study of 
linkages with transverse coupler link vibrations which are uncoupled with 
the vibrations of the rest of the system. They recommend relatively robust 
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mechanisms for high speed linkages. However, robust mechanisms carry 
weight penalties which are unacceptable in many applications. As demand 
for faster systems is increasing and because vibration is the phenomenon 
which most often limits the speed at which a high-speed mechanism can 
operate, it is important to develop tools for predicting frequency response 
for mechanisms consisting of flexible linkages and elastic input and output 
shafts. To date only a few authors have included in their models the effects 
of flexibility in the system input and output shafts or the effects of inertial 
loads on the ends of these shafts. In 1967, Capellen [13] considered input 
and output shafts but did not account for any elasticity in the linkage it- 
self. In 1977, Kohli, Hunter and Sandor [50] included the elastic effects of 
supports and input and output shafts on the vibration of an elastic linkage; 
they applied fourier series representations to model the deflections of the 
links. In this paper we consider the case in which a four-bar linkage with an 
elastic rocker or output link is driven by a motor through a flexible input 
shaft and drives, through connection to the rocker link, a flexible output 
shaft which in turn drives an inertial load. (Figure C.l). 

Starting with a continuous parameter model of the system, the behavior 
of the rocker-beam is analyzed, and from these results a lumped parame- 
ter model is constructed. Then the equations of motion are linearized to 
get a system of differential equations with periodic coefficients. Finally a 
dimensional analysis leads to a procedure for design of such linkages whose 
resonances do not match frequencies of system energy input. 
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C.3 Description of the Continuous Parame- 
ter Model 

Figure C.l shows a crank-rocker mechanism with flexible input and out- 
put shafts. Two discrete and two continuous generalized coordinates are 
required to describe this system. They are the following: 

• 0o , the angular position of the motor rotor 

• 9\ , the angular position of the input link 

• u(y, t) , the displacement from static equilibrium of all points on the 
rocker 

• <f>(x,t) , the torsional displacement from static equilibrium of all 
points on the output shaft. 

We shall derive the equations of motion of this continuous system for dis- 
crete configurations; that is, for small motions about a given position of 
the crank; using the energy Method [20]. 

C.3.1 Kinetic co-energies 

The system elements have the following kinetic co-energies: 
. Output shaft: 7? = \ f L > Pa I p .(%) 2 dx 
. Inertia / 4 : T; = i/ 4 (f )£. 
. Rocker: T^ = \ ^ p z A,{^dy 

• Crank: T 4 * = \hO\ 

• Motor inertia: T 5 * = A/,$ 



APPENDIX C. PRELIMINARY STUDY OF VIBRATIONS 203 

• Coupler: T * = ^AO*, where 

A — l2, v i\ + m 2 (Lj + -£i\ + LiL 2 t2 cos (#i — #2)) where t 2 and 13 (used 
later) are the angular velocity ratios as defined by Paul [60]. 

V _ 02. = Li »in(fl 3 - g|) 
12 0, L 7 sin (0-J (?;,) 

,• _ ja. _ Li gia(ga-Oi) 
3 By t 3 8in(flj-fl 3 ) 

• Total kinetic energy: T* = T; + T 2 + ... + T 6 * 

C.3.2 Potential Energies 

The system elements have the following potential energies: 

• Output shaft: Vk = £ / i3 <?./,. (§£)*dx 
. Rocker: V 2 = j /*> ^/»(g|)*rfy 

• Input shaft: V 3 = §#i(0i - O ) 2 

• Coupler axial spring: V 4 = \K 2 (NiLn0i — iV 2 Ui) 2 , where 

Ni — sin (#i — 6 2 ) and N 2 = sin (#3 — $2) are the factors that project 
the displacements at the ends of the coupler in the direction of the 
coupler, and 

ui = u(L, t). 

• Total potential energy : V = Vl + V2 + V3 + V4 

C.3.3 Generalized Forces 

Variation of virtual work is 6W = T 66 - CxK^Oi - C 2 (^) x=0 6iJ>o. 
Then Z 0o = fj^ = T , Z 6l = -Crfi, Z^ = -C 2 (-^) x=0 
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Using Hamilton's principle [20] we obtain: 

VI = (6L + J2 ZjS^dt, L = T* -V 
Jtl j=i 

Substitution of above expressions for T and V into the expression for the 

Variational indicator and then integration by parts both timewise and 

spacewise produces the equations of motion and boundary conditions at 

the vanishing of the Variational indicator for arbitrary variations of the 

generalized coordinates. These are as follows: 

P*& + G/]£=0in0<x<L, 
Ps^gr + £ 3 / 33 gf = OinO < y < L z 
h i ££+GJ p &=0atz = L. 

GJ,M X = ~ <*SLo + *W»(&).-o = Oaf x = 0,y = 

p = 0aty = L 3 

(Ji + A)h + C& + Ktfa - O ) + KiNiL^NyLtfx - N 2 u L ) = 0, 

at 0i 

J £o + ^(00-01) = Oaf O 

£3/33 g? +K a N i {N 1 L 1 e l - N 2 u L ) = 0aty = L z 

For this derivation the equation of motion for 0i has been linearized 
by considering small motions about a given configuration. To solve the 
eigenvalue problem, we put temporarily C\ = C 2 = 0, Tq = and try 
solutions of the form: 

ip(x, t) - a n (x) sin (u„t + <p) 
"(y. t) = b n (y) sin (w„f + <f>) 
9i(t) = 0x sin (u n t + 4>) 
tfo(0 = ^osin(w„« + ^) 
The substitution of these solutions into the above equations gives: 

p„uila u {x) + G,—^ = 0,in0 < x < l 3 (C.80) 

(71" 
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- p 3 A 3 ulb n {x) + £3/33 34 = 0, inO < y < L 3 {CM) 

-ul{h+A)e i + K i {0 i -h)+K 2 N*L\0 i -K 2 N i N 2 L l b n {L)=0,at0 l (C.82) 

- I Q u% + K x (0 O - 9i) = 0, at O (C.83) 

GJ P , -j± - I 4 u 2 n a n = 0,atx = L s {CM) 

GJ P . ^ + £3/33 ^^ = 0, at x = 0, y = (C.85) 

5 2 6„ft/) 

-^P = 0,a*y = L 3 (C.86) 

q3i 

- £3/33 yf + K 2 Nlb n {y) - K^N^LA = 0,aty = L 5 (C.87) 

Also from geometry: 6 n (0) = 0, and a n (0) = (^ a ) tf =o 
From (1) we get: 

a n (i) = Ci sin —= + C 2 cos -= 

From (2) we get: 

K{y) = £3 sin (3y + C 4 cos /?y + C 5 sinh /?y + C cosh /?y 

where /? = (>-%gY>* 

After some manipulations in order to eliminate the constants C\ through 
Co, 0i and 0o we can get the transcendental equation: 

Z x = Z 2 (CSS) 

where 



Z 1 = 



ff^Uanfl&a + 2sinh/3£ 3 ' 



„ , /F7T G.Io,+/G./p,coe( w " f '»- )-/ 4 u) n 8in( ""A', ) 

0{cot0L 3 - cot/i0Z 3 ) + 2^ fi^V G - /p ' [ V , V°- f "- l^p- 

V ; K G ' J " w - l G./ (> . v /G./^a n (-^fefc) + /4-„co.(^^) 

and 

Z = S ^ K ^ - ^^ ~ Iffl) +P* l { (cot0L 3 -coth/?L 3 ) 
2 -/? 3 (cot/?L 3 cos5L 3 + sm/?L 3 + coth /?L 3 cosh /3Z, 3 -sinh/?L 3 ) 
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with 

-K X + \-»l{h + A) + K i + KiNfLWi*^) 
K^NiNiL). 
From equation C.88 we can obtain an infinite set of eigenvalues which con- 
sists of a subset reflecting rocker bending, a subset reflecting shaft torsional 
vibration and a subset reflecting lumped parameter vibration modes. The 
behavior of the rocker is particularly important because one of its eigenval- 
ues may be at or near the system energy input frequency for some system 
configurations and the magnitudes of its eigenvalues depend on all the pa- 
rameters of the system. For example, if the stiffness of the input shaft is 
large and Iq (motor inertia) is large, the rocker tends to behave as a beam 
whose upper end is pinned, whereas if the input shaft stiffness is small, the 
rocker tends to behave as a free-end beam. Additionally, if both the stiff- 
ness of the output shaft and the output inertia I 4 are large, then near the 
fixed bearing the rocker tends to behave aa a clamped-end beam while if 
this stiffness is small, it behaves as a pinned-end beam. To illustrate this 
behavior we define a parameter o as follows: 



2 / ^3^33 
<*>!, tending = a 



p 3 A z Lj 

where a depends on the boundary conditions. Values of a for typical bound- 
ary conditions are as follows [29]: 

• hinged-clamped a = 3.91 

• hinged-hinged a = ir 

• cantilever a = 1.87 

• free-hinged a = 0. 

After solving the eigenvalue problem for several values of Ki and K , for 
fixed K s we obtain the behavior of a vs Ky for the two values of K s (sec 
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Figure C.2: Behavior of a vs. K\ for various values of K, 

Figure C.2). Note that in the case where both input and output shafts are 
soft, a for the rocker takes a value (zero) appropriate for a beam hinged or 
pinned at its base and free at the coupler end. When both shafts are stiff, 
the value of a is that of a beam clamped at the base and hinged at the 
other end. For a soft input and stiff output shaft, the value of a is that of 
a cantilever , while for a stiff input and a soft output shaft a has the value 
fitting a hinged-hinged beam. 



C.4 Description of Simplified Models 

For most actual systems, the motor inertia is very high compared to the 
inertias of each link of the mechanism and the input shaft is short and 
therefore very stiff. Existence of these conditions implies that the rocker 
will behave near the rocker-coupler joint as a pinned beam. Based on this 
implication, a bending deformation pattern for the rocker can be assumed 
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Figure C.3: Lumped parameter model of the system 

and a much simpler system model, described by a set of ordinary differential 
equations, can then be developed. Validity of this assumption , for a specific 
system can be checked by comparing results of this simplified analysis to 
those obtained from the analysis of the more precise continuous model. 
From the simple model, a finite number of eigenvalues and mode shapes 
will be obtained.For most problems, only a finite range of eigenvalues is 
of interest and so the loss of higher eigenvalues is not a serious limitation, 
and the simplified system equations are much easier to solve than are those 
for the continuous model. For the case of a large motor inertia and a stiff 
input shaft, the vibration induced deflection pattern of the rocker link can 
be assumed to be sinusoidal. 



«(y. t) = <j{t) sm — 

For this case, a nondimensional generalized coordinate z = q/Lz is chosen. 
The output shaft is modeled as a masslcsa spring K a . Figure C.3 shows 
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this simplified model. The Lagrangian for this system is as follows: 

L = \6][h + I^ 9 i\ + m 2 {L\ + %i\ + L 1 L 2 i 2 cos {0, - 2 ))] + 
±/„0 o 2 + {hb\ + ±m 3 (f i\9\ + %z 2 + ^ijxz) [CM) 

-J#i(*i - O ) 2 - \K.{0 A -9 3 - *zf - "i^z 2 
where wz, the slope of the rocker beam near the fixed bearing relative 
to its position at static equilibrium, is considered approximately equal to 
the angle made by a tangent to the beam at the fixed bearing with respect 
to the position of the beam in static equilibrium. From the Lagrangian 
we can derive equations of motion which are nonlinear ordinary differential 
equations. Note that the term in 7713 is the kinetic energy of the beam and 
is obtained by 

KE = -p 3 A 3 I v 2 dy 
. z Jo 

, where 

Try 

v = 9 3 y + zL 3 sin — 

L 3 

and the term in z is the potential energy of the beam and is obtained by 

To obtain a more precise but still a simple set of system equations, we 
can, instead of taking zL 3 aia(^-) as the first natural mode of the beam, 
take the actual mode shape obtained from the continuous parameter model 
discussed above, then 

&i (y, = zL 3 [C 3 sin (3y + C 4 cos /3y + C 5 sinh /?y + C 6 cosh /3y] 

where C 3 ,C t ,Cs and C& are determined for a particular case from the con- 
tinuous model. The kinetic and potential energies then become: 

KE = \m 3 {&'el + hz* + f 2 b 3 z) 

V = \^fzz 2 
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Figure C.4: Variation of three natural frequencies with crank position for 
quasi-static case 

where /i, ft, /j are constants derived for the particular case. When these 
energies are incorporated into the Lagrangian (10), we obtain equations of 
motion, this time by use of Lagrange's equations, which more precisely 
describe the bending motion of the link, (see equation (11) which for w = 
gives the quasi-static case.) For a given linkage configuration we linearize 
the equations of motion and since the coefficients are constant we can form 
and solve the eigenvalue problem. When this solution is obtained for various 
configurations we can determine the variation of the system eigenvalues in 
one cycle of the crank as a function of crank angle. Figure C.4 shows 
how the eigenvalues change for discrete crank positions. Shown are results 
for the continuous model. Results obtained from the lumped mass models 
are very similar. To obtain values for these eigenvalues, we have used the 
system parameters calculated from drawings of a linkage which is part of 
an industrial sewing machine that has presented its manufacturer with a 
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serious vibration problem during high speed operation. 

This system exhibited a strong vibration at 400 Hz with the crank speed 
at 100 Hz, illustrating that even for a small ratio of crank length to ground 
link length, 0.13 in this case, the geometric nonlinearity brought about an 
important parametric excitation. Values of the chosen system parameters 
are shown in table C.l. 

The lumped parameter system has four generalized coordinates and has 
therefore four eigenvalues, one of which is zero. Note that for the frequency 
range of interest (0-1 000 Hz) the eigenvalues are nearly constant around the 
cycle and this allows us to consider fuf^th to be constants. 

C.5 Comparison of the Continuous and Lumped 
Mass Models 

In order to illustrate the comparison of the eigenvalues obtained from the 
lumped parameter models with the ones obtained from the continuous pa- 
rameter model which are known to be more precise, we have solved each 
eigenvalue problem for various system parameters with the crank position 
fixed at 90 degrees . These results are shown in figure C.5 and Figure C.6. 
Differences in the first eigenvalues are seen to be small for a wide range of 
input shaft stiffnesses for two choices of output shaft stiffnesses. Maximum 
difference is 5 % of the peak value. 

C.6 Linearized System Vibration for Steady- 
Motor Rotation 

In order to linearize the equations of motion we assxime that the motor 
rotates at constant speed and that a small amplitude vibration can be 
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L { = .02m 


mx = 4.08341£ - 2 Kg 


L 2 = .035m 


m 2 = 2.G9775£ - 2 Kg 


L 3 = -150m 


m 3 = 1.70105£ - 2 Kg 


G = .155m 


To = 2.82528JS - lKg - m 2 


7 33 = 9.037£ - 12m 4 


h = 5.12082^ -6Kg-m 2 


E z = 2.0601^1 LN/m 2 


h = 2.74C8J57 - 6iiCg - m 2 


Ki = 5258 N-m/rad 


7 3 = 1.2753/5 - 4Kg - m 2 


K„ = 202 N-m/rad 


7 4 = 1.56£ - 3Kg - m 2 



Table C.l: Data of a Particular Mechanism 
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Figure C.5: Comparison of behavior of first natural frequency in bending 
obtained from the continuous and the lumped models. Case of jj:*- = 0.1 
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Figure C.6: Comparison of behavior of first natural frequency in bending 
obtained from the continuous and the lumped models. Case of jp- = 10. 

superimposed on the rigid body motion. Then we can obtain a system of 
linear differential equations with periodic coefficients. To do this we define 
a new generalized coordinate V>i which will replace O and 9 X ,thus reducing 
the number of degrees of freedom to three, as follows: 

ip l =8 l - O , 0i = V>i + w* 

where u; is the angular velocity of the motor and t is the time. To first 
order approximation 

03 = #3 + *S^1 

where 03 is the rocker position for rigid body motion and 3 is the corrected 
value. Also 

0! = V>! + W, 

0! = tf lt 

3 = i 3 U> + t 3 ^! 
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The equations of motion obtained from this system can be linearized as 

follows: 

V>i = Stpi (already small) 

z = Sz (already small) 

#4 = &3 + S64 

These equations then are : 



A £m 3 L§i 3 

I 4 

£m 3 i#3 ° /l^3^3 



> + 



Stpi 
S9 4 
Sz 

TC2H 
-TC 2 



[A'u + Ci. - C 2 i\ -C 2 i z 
—C 2 i z C 2 

TC 2 i 3 -TC 2 T 2 C 2 + C Z 

K\ + K 2 i z —K 2 i$ TK 2 %3 

—K 2 i z K 2 —TK 2 

TK 2 i z -TK 2 T 2 K 2 + f 3 2g* 

^ W 2 -C 2 t> 

— /4.K3W 2 + C 2 l3(jJ 

^■m z L\K' z u 2 - C 2 i z w 



60 4 
Sz 
Sip! 

se A } = 

Sz 



> + 



1 



where 



(C.90) 



• A = h + h C9 i\ + m 2 {L\ + \L\i\ + L x L 2 i 2 cos (B x - 9 2 )\ 

• T = a(C z + C$), obtained from mode shape(continuous model) 

• #3 = k ii? 

» u at 

These equations (11) can be integrated numerically to get the steady state 
response and determine the amplitudes of vibration and the frequency con- 
tent. 

Because the process of numerical integration can be expensive in terms 
of computer time, it is helpful to perform a dimensional analysis before 
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integrating. Then each numerical solution will apply to a family of mecha- 
nisms. Define the nondimensional parameters: 



E3 /33 



D 2 = ^& 



D, = c '^ 



D* 



_Ca_ 



2y/K.u ' ~ 3 2v//C 3 «3 

and the following parameters with dimensions of frequency: 



n 
n 4 



Ql ~ ys£tf» n 3 - y sfei 



n * = yJmHu*' n 8 = y s5t»» n »o = y dSj 

Define also a dimensionless time t = n 4 t Then, 

The equations of motion (11) become after some manipulations: 




A 




(&) 2 







2*3 


6ti 







h 



+ 



2*3 

&& + 2D & k ~ ^(ft) 1 ] -2*3^2(8t) 2 

-2*3iM{£) 2 

2Tt 3 Z? 2 (g) 2 

(§})+«•.(&)■ -*3(8^) 2 



+ 



2Ti 3 D 2 (%)> 

2D 2 (£) 2 -2rD 2 (g) 2 

-2z- 3J D 2 (^) 2 2r 2 A(^) 2 ^+ %D 3 %? 



<S0 4 

<5z 



® 2 

£ (£) 2 - 2t!JM&)'(£) 



r.-,(jfe)» 






<J0 4 
6* 






-(etm(nT) 2 + 2i 3J D2(et) 2 (nT) 
[ ^* 3 (£) 2 - 2^ 2 (ft) 2 (^) 



(C.91) 



where 
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Figure C.7: Frequency response of two members of a family of mechanisms. 
The original data corresponds to ordinate 2 
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3 



• * = (ig) 2 + t(%) 2 + (£) 2 + (£) 2 < 2 cos (e l -e 2 ) + 'i 

• T = a(C' 3 + C5), obtained from mode shape (continuous model) 

<> u at 

From these equations (12) we can observe clearly the consequences of dy- 
namic similarity [64]. If for example the system has a small value for the 
damping ratio as most mechanisms do, then changing all the stiffnesses 
proportionally, or all the inertias proportionally will not change the dy- 
namic response expressed in nondimensional variables. Therefore, all the 
systems that have the same ratio of j£, jjSj*,^, ^f-> j* belong to a fam- 
ily whose dynamic response is unique. By translating the results from the 
nondimensional space to a dimensional space, we can produce the following 
graph. (See figure C.7). 
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This figure shows the frequency response for two members of the family. 
Each data point represents an angular displacement of the inertial mass on 
the output shaft as determined from a numerical integration with respect 
to its expected position for rigid body motion during operation at the given 
motor angular velocity. The angular displacement is an RMS displacement 
measured at 256 points during six cycles of the input crank after the system 
had reached equilibrium. For the two members of the family shown, all 
system stiffnesses differ by a factor of 2; damping ratios and inertias for both 
are the same. From the figure we see that if amplitudes are plotted vs w 2 , 
points of constant response amplitude for the two family members occur at 
frequencies which lie on a straight line emanating from the origin. This is as 
expected on the basis of the system equations (12) and of dynamic similarity 
[64]. After the frequency response for one member of the family is obtained, 
these results can be extended to any other member simply by scaling the 
frequency axis appropriately. Having performed this analysis for a family 
from which a designer expects to pick a member, and having determined 
from design considerations the frequencies (driving frequency for example) 
at which energy will be introduced into the system, the designer can now 
pick a member of the family whose frequency response at the known input 
frequencies is low. 

C.7 Summary 

Analysis of a four-bar linkage system with flexible input and output shafts 
has shown that for typical system parameters, the flexibility in the driv- 
ing and driven shafts have a significant influence on the system response. 
Quasi-static response has been reported for a continuous parameter system 
model, and this response has been compared to that obtained for much 
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simpler lumped parameter system models. For wide ranges of system pa- 
rameters, the lumped system responses are similar to those obtained from 
the more precise continuous system model. Values for first and second 
system natural frequencies are shown to vary only slightly with changes in 
system geometry. Two sets of nondimensionalized lumped parameter model 
system equations, one assuming sinusoidal mode shape for the rocker link 
and one using for rocker link mode shapes the shapes determined from 
the quasi-static analysis of the continuous model, have been integrated nu- 
merically to establish dynamic system response. Resulting plots of RMS 
displacement of the output shaft about its calculated rigid body motion po- 
sition vs. frequency of the input shaft rotation can be examined, and from 
these results, and on the basis of dynamic similarity, design parameters for 
the system can be chosen so as to minimize RMS response at frequencies 
to which the system is likely to be subjected. 
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